The Museum of HP Calculators


Miscellaneous Mathematical Functions for the HP-41

This new version updated 12/12/05. New functions and Improvements noted in red. Previously updated 12/13/04.

This program is Copyright © 2004 by Jean-Marc Baillard and is used here by permission.

This program is supplied without representation or warranty of any kind. Jean-Marc Baillard and The Museum of HP Calculators therefore assume no responsibility and shall have no liability, consequential or otherwise, of any kind arising from the use of this program material or any part thereof.


Overview
 

  1°)   Gamma , Digamma and Polygamma Functions ( real arguments ) ; Gamma Function in the complex plane ; Log Gamma(x) for large arguments.1/Gam(x)
          and Digamma Function in the complex plane ( new )
  2°)   Bessel & modified Bessel Functions of the first and second kind. Integral of Jn(x)
  3°)   Kummer's ( or Confluent Hypergeometric ) Functions  M(a;b;x)
          and  M(a;b;z)  i-e complex variable  ( new )
  4°)   Parabolic Cylinder Functions  U(a;x)  V(a;x)  W(a;x)
  5°)   Hypergeometric Functions  F(a;b;c;x)   ( 2 programs )
  6°)   Legendre Functions ( Legendre Polynomials & Associated Legendre Functions of the first and second kind )  ( 7 programs in all )
  7°)   Laguerre Polynomials
  8°)   Hermite Polynomials
  9°)   Chebyshev Polynomials ( of the first and second kind ) & the "Connaissance des Temps"
 10°)  Other Orthogonal Polynomials:  ( new )
          Ultraspherical , Generalized Laguerre , Jacobi
 11°)  Struve Functions
 12°)  Theta Functions
 13°)  Riemann Zeta Function
 14°)  Weierstrass Elliptic Functions ( real arguments & complex arguments )
 15°)  Exponential, Sine and Cosine Integrals  ( Ei(x) , Si(x) , Ci(x) , Shi(x) , Chi(x) )  ( slightly improved )
          and  En(x) ( new )
 16°)  Fresnel Integrals
 17°)  Error Function  ( slightly improved )
 18°)  Coulomb Wave Functions  ( new )
 19°)  Debye Functions  ( new )
 20°)  Dawson's Integrals  ( new )
 

-Some of these functions are obtained by a series expansion and, sometimes, it will yield a good accuracy only if x is not too "large"
  but it may also depend on the other parameters.
-Other programs use asymptotic expansions ( for large arguments ). Though usually divergent, these sums give an error smaller than the first omitted term.
 

1°) Gamma & Polygamma Functions
 

        a) Gamma Function

-The asymptotic formula  e-x xx-1/2 (2pi)1/2 ( 1 + 1/(12x) -1/(360x3) ) is used if x > 16
-Otherwise,  an integer n is added to x such that  x+n > 16  and this formula is then applied to x+n
  with the relation:  Gam(x+n) = (x+n-1) (x+n-2) ...... (x+1) x Gam(x)
-Gam(x) is defined if  x # 0 ; -1 ; -2 ; .......
-Synthetic register M may be replaced by any unused register.
-This short program is called as a subroutine by several other programs hereafter.
 

Data Registers:  /
Flags: /
Subroutines: /
 

01  LBL "GAM"
02  FRC
03  X#0?
04  GTO 00
05  LASTX
06  1
07  -
08  FACT
09  RTN
10  LBL 00
11  16
12  LASTX
13  STO M
14  GTO 02
15  LBL 01
16  1
17  +
18  ST* M
19  LBL 02
20  X<Y?
21  GTO 01
22  ENTER^
23  X^2
24  SIGN
25  LASTX
26  30
27  *
28  1/X
29  -
30  X<>Y
31  12
32  *
33  /
34  E^X
35  0
36  X<> M
37  /
38  X<>Y
39  1
40  E^X
41  /
42  R^
43  Y^X
44  *
45  X<>Y
46  PI
47  *
48  ST+ X
49  SQRT
50  *
51  END

( 69 bytes / SIZE 000 )
 
 
      STACK         INPUT      OUTPUT
           X              x        Gam(x)

 

Example:    Compute   Gam(PI)  and  Gam (-6.14)

      PI   XEQ "GAM"  yields     2.288037783    ( in 5 seconds )
  -6.14      R/S            ------    -0.007872567    ( in 7 seconds )
 

        b) Digamma Function
 

-The Digamma ( or Psi ) function is defined by   Psi(x) = Gam'(x)/Gam(x)   where  Gam'(x) is the first derivative of Gam(x)
-The asymptotic expansion  Psi(x) = ln x - 1/(2x) -1/(12x2) + 1/(120x4) - 1/(252x6) + 1/(240x8)  is used for  x > 8
  together with the property:  Psi(x+1) = Psi(x) + 1/x
 

Data Registers:  /
Flags: /
Subroutines: /

-Synthetic registers M & N may be replaced by R00 & R01
-In this case, delete line 41 and replace line 02 with    0   STO 00   RDN
 

01  LBL "PSI"
02  CLA
03  8
04  X<>Y
05  LBL 01
06  1/X
07  ST+ M
08  X<> L
09  1
10  +
11  X<Y?
12  GTO 01
13  1/X
14  STO N
15  X^2
16  ENTER^
17  STO Z
18  20
19  /
20  21
21  1/X
22  -
23  *
24  .1
25  +
26  *
27  1
28  -
29  *
30  12
31  /
32  RCL N
33  LN
34  LASTX
35  2
36  /
37  +
38  -
39  RCL M
40  -
41  CLA
42  END

( 61 bytes / SIZE 000 )
 
 
 
      STACK         INPUT      OUTPUT
           X              x        Psi(x)

Example:    Compute Digamma(-1.6) & Digamma(1)

  -1.6  XEQ "PSI"  >>>  Digam(-1.6) = Psi(-1.6) = -0.269717876   ( in 4 seconds )
     1         R/S        >>>   Psi(1) = -0.577215665 = the opposite of the Euler's constant  ( 0.5772156649 )
 

        c) Polygamma Functions
 

-The following program calculates the nth derivative of the Psi function   ( n = 0 ; 1 ; 2 ; ..... )
-  Psi'(x) = the Trigamma function ,  Psi''(x) = the Tetragamma function ... etc ...
-The asymptotic expansion of the Psi-function is derived n times and the recurrence relation  Psi(n)(x+1) = Psi(n)(x) + (-1)nn! x-n-1 is used.
 

Data Registers:  /
Flags: /
Subroutines: /

-Synthetic registers M N O may be replaced by R00 R01 R02
-In this case, delete line 95 and replace line 02 with    0   STO 00   RDN

01  LBL "PSIN"
02  CLA
03  X<>Y
04  STO O
05  8
06  +
07  X<>Y
08  STO N
09  LBL 01
10  RCL O
11  1
12  ST+ N
13  +
14  CHS
15  Y^X
16  ST+ M
17  CLX
18  RCL N
19  X<Y?
20  GTO 01
21  1/X
22  X^2
23  STO Y
24  RCL O
25  7
26  +
27  FACT
28  7
29  FACT
30  /
31  *
32  20
33  /
34  RCL O
35  5
36  +
37  FACT
38  2520
39  /
40  -
41  *
42  RCL O
43  3
44  +
45  FACT
46  60
47  /
48  +
49  *
50  RCL O
51  1
52  +
53  FACT
54  -
55  *
56  12
57  /
58  RCL N
59  RCL O
60  Y^X
61  ST/ Y
62  RCL N
63  *
64  RCL O
65  FACT
66  ST* M
67  X<>Y
68  ST+ X
69  /
70  -
71  RCL O
72  X=0?
73  GTO 02
74  1
75  -
76  FACT
77  RCL N
78  RCL O
79  Y^X
80  /
81  -
82  GTO 03
83  LBL 02
84  X<> N
85  LN
86  +
87  LBL 03
88  RCL M
89  -
90  1
91  CHS
92  RCL O
93  Y^X
94  *
95  CLA
96  END

( 138 bytes / SIZE 000 )
 
 
 
      STACK        INPUTS     OUTPUTS
           Y             n            /
           X             x       Psi(n)(x)

                                   ( if n = 0  we have the Digamma ( or Psi- ) function
                                     if n = 1  ------------ Trigamma function
                                     if n = 2  ------------ Tetragamma function and so on ... )

Example:    Calculate Digamma(-1.6)  Trigamma(-1.6)  Tetragamma(-1.6)  Pentagamma(-1.6)

   0   ENTER^
-1.6  XEQ "PSIN"  >>>  Digam(-1.6) = Psi(-1.6) = -0.269717877

   1   ENTER^
-1.6     R/S             >>>  Trigam(-1.6) = 10.44375936

   2   ENTER^
-1.6     R/S             >>>  Tetragam(-1.6) = -22.49158811

   3   ENTER^
-1.6     R/S             >>>  Pentagam(-1.6) = 283.4070827    ... etc ...   ( execution time is about 13 seconds for these examples )
 

        d) Gamma Function in the complex plane

-This program employs the same asymptotic formula as "GAM"

Data Registers:  R00 to R05: temp
Flags: /
Subroutine: /
 

01  LBL "GAMZ"
02  RAD
03  STO 01
04  X<>Y
05  STO 02
06  1
07  STO 03
08  CLX
09  STO 04
10  LBL 01
11  16
12  RCL 01
13  X>Y?
14  GTO 02
15  RCL 02
16  X<>Y
17  R-P
18  ST* 03
19  X<>Y
20  ST+ 04
21  1
22  ST+ 01
23  GTO 01
24  LBL 02
25  RCL 02
26  RCL 01
27  R-P
28  2
29  CHS
30  ST* Z
31  Y^X
32  360
33  CHS
34  /
35  P-R
36  12
37  1/X
38  +
39  R-P
40  RCL 02
41  RCL 01
42  R-P
43  RDN
44  ST- Z
45  X<> T
46  /
47  P-R
48  X<>Y
49  RCL 02
50  -
51  STO 00
52  X<>Y
53  RCL 01
54  -
55  STO 05
56  RCL 02
57  RCL 01
58  .5
59  -
60  R-P
61  RCL 02
62  RCL 01
63  R-P
64  LN
65  R-P
66  RDN
67  ST+ Z
68  X<> T
69  *
70  P-R
71  X<>Y
72  RCL 00
73  +
74  RCL 04
75  -
76  X<>Y
77  RCL 05
78  +
79  E^X
80  RCL 03
81  /
82  PI
83  ST+ X
84  SQRT
85  *
86  P-R
87  DEG
88  END

( 113 bytes / SIZE 006 )
 
 
      STACK        INPUTS      OUTPUTS
           Y             y             B
           X             x             A

  where  Gam(x+iy) = A + iB

Example:    Compute   Gam(1+2i)

     2   ENTER^
     1   XEQ "GAMZ"   yields  0.151904003
              X<>Y                       0.019804881        thus,   Gam(1+2i) = 0.151904003 + 0.019804881. i      ( in 21 seconds )

Note:   -If you prefer to obtain  Ln(Gam(x+iy)) , simply replace line 86 by LN
           -For instance,  Ln(Gam(2+3i)) = -1.876078781 + 0.129646321 i
 
 

       e) Logarithm of the Gamma Function ( large arguments )
 

-This program evaluates  log | Gam(x) |  and it's especially useful if n > 69

Data Registers: /
Flags: /
Subroutines: /
 
 

01  LBL "LOGAM"
02  ENTER^
03  ABS
04  LOG
05  STO M
06  16
07  RCL Z
08  LBL 01
09  1
10  +
11  STO Z
12  ABS
13  LOG
14  ST+ M
15  X<> Z
16  X<Y?
17  GTO 01
18  ENTER^
19  X^2
20  SIGN
21  LASTX
22  30
23  *
24  1/X
25  -
26  X<>Y
27  12
28  *
29  /
30  X<>Y
31  -
32  10
33  LN
34  /
35  X<>Y
36  ENTER^
37  LOG
38  *
39  +
40  X<>Y
41  PI
42  *
43  ST+ X
44  SQRT
45  LOG
46  +
47  0
48  X<> M
49  -
50  END

( 72 bytes / SIZE 000 )
 
 
      STACK        INPUTS      OUTPUTS
           X              x     log | gam(x) |

Example:    1000  XEQ "LOGAM"  yields   log | gam(1000) |  = 2564.604644   whence  gam(1000) = 4.02387  102564
 

-If you prefer   ln | gam(x) |  instead of  log | gam(x) |  , delete lines 32-33-34 and replace all the LOG  with  LN
 

       f) Reciprocal of the Gamma Function
 

-It's sometimes worthwhile to use a program which computes 1/Gam(x) because this function is defined for any x-value.
-Thus, it may avoid tests when  x = 0 ; -1 ; -2 ; -3 ; ......   ( cf  6°) g) )

Data Registers: /
Flags: /
Subroutines: /
 

01  LBL "1/G"
02  STO M
03  16
04  X<>Y
05  GTO 02
06  LBL 01
07  1
08  +
09  ST* M
10  LBL 02
11  X<Y?
12  GTO 01
13  ENTER^
14  ENTER^
15  X^2
16  30
17  *
18  1/X
19  1
20  -
21  X<>Y
22  12
23  *
24  /
25  E^X
26  0
27  X<> M
28  *
29  X<>Y
30  1
31  E^X
32  /
33  R^
34  Y^X
35  /
36  X<>Y
37  PI
38  *
39  ST+ X
40  SQRT
41  /
42  END

( 59 bytes / SIZE 000 )
 
 
 
      STACK        INPUTS      OUTPUTS
           X              x        1 / gam(x)

Example:      PI  XEQ "1/G"  >>>>  1/Gam(pi)  =  0.437055720  ( in 5 seconds )
                     -3       R/S         >>>>   1/Gam(-3) = 0
 
 

        g) Digamma Function in the complex plane
 

Data Registers:  R00 thru R05: temp
Flags: /
Subroutines: /

01  LBL "PSIZ"
02  STO 01
03  X<>Y
04  STO 02
05  CLX
06  STO 03
07  STO 04
08  LBL 01
09  10
10  RCL 01
11  X>Y?
12  GTO 02
13  RCL 02
14  X<>Y
15  R-P
16  1/X
17  X<>Y
18  CHS
19  X<>Y
20  P-R
21  ST+ 03
22  X<>Y
23  ST+ 04
24  1
25  ST+ 01
26  GTO 01
27  LBL 02
28  RCL 02
29  RCL 01
30  R-P
31  X^2
32  STO 00
33  1/X
34  X<>Y
35  ST+ X
36  STO 05
37  CHS
38  X<>Y
39  21
40  /
41  P-R
42  .1
43  -
44  R-P
45  RCL 00
46  /
47  X<>Y
48  RCL 05
49  -
50  X<>Y
51  P-R
52  1
53  +
54  R-P
55  12
56  /
57  RCL 02
58  RCL 01
59  R-P
60  RDN
61  ST- Z
62  X<> T
63  /
64  P-R
65  .5
66  +
67  R-P
68  RCL 02
69  RCL 01
70  R-P
71  RDN
72  ST- Z
73  X<> T
74  /
75  P-R
76  RCL 02
77  RCL 01
78  RAD
79  R-P
80  LN
81  DEG
82  R^
83  ST- Z
84  X<> T
85  -
86  RCL 04
87  ST- Z
88  X<> 03
89  -
90  END

( 118 bytes / SIZE 006 )
 
 
 
      STACK        INPUTS      OUTPUTS
           Y             y             b
           X             x             a

 with  Psi(x+i.y) = a + i.b

Example:     0.7  ENTER^
                    1.6  XEQ "PSIZ"  >>>>  0.276737830    X<>Y    0.546421305    ( in 23 seconds )

  Whence   Psi(1.6+0.7.i) = 0.276737830 + i. 0.546421305
 

2°) Bessel Functions
 

        a) Bessel & modified Bessel Functions of the first kind  Jn(x) &  In(x)

Formulae:       Jn(x) = (x/2)n [ 1/Gam(n+1)  +  (-x2/4)1/ (1! Gam(n+2) )  + .... + (-x2/4)k/ (k! Gam(n+k+1) ) + ....  ]      n # -1 ; -2 ; -3 ; ....

                      In(x) = (x/2)n [ 1/Gam(n+1)  +  (x2/4)1/ (1! Gam(n+2) )  + .... + (x2/4)k/ (k! Gam(n+k+1) ) + ....  ]      n # -1 ; -2 ; -3 ; ....

Data Registers:  /
Flag:  F01
Subroutine:  "GAM"  ( Gamma Function )

-Synthetic registers N & O may be replaced by any unused data registers.
 

01  LBL "JNX"
02  2
03  /
04  STO N
05  X<>Y
06  STO O
07  1
08  ENTER^
09  STO M
10  ENTER^
11  LBL 01
12  X<> Z
13  FC? 01
14  CHS
15  RCL N
16  X^2
17  *
18  RCL M
19  ST/ Y
20  RCL O
21  +
22  /
23  ISG M
24  CLX
25  ST+ Y
26  X<> Z
27  X#Y?
28  GTO 01
29  RCL N
30  RCL O
31  Y^X
32  *
33  X<> O
34  1
35  +
36  XEQ "GAM"
37  RCL O
38  X<>Y
39  /
40  CLA
41  END

( 70 bytes / SIZE 000 )
 
 
 
      STACK        INPUTS     OUTPUTS
           Y             n        f(x)
           X             x        f(x)

 where  f(x) = Jn(x)  if flag F01 is clear  and  f(x) = In(x)  if flag F01 is set.

Example:    Calculate   J0.7(1.9)  and  I0.7(1.9)

   CF 01
   0.7   ENTER^
   1.9   XEQ "JNX"   yields   J0.7(1.9) = 0.584978102   (  in 12 seconds )

   SF 01
   0.7   ENTER^
   1.9   XEQ "JNX"   yields   I0.7(1.9) = 1.727630598   (  in 12 seconds )
 

Notes:  if  n  is a positive integer,  J-n(x) = (-1)n Jn(x)   and    I-n(x) = In(x)
         J0(0) = 1 but this program doesn't work if x = n = 0 ( DATA ERROR line 31 )
    if you add  X=Y?  X#0?  GTO 00  SIGN  X<>Y  LBL 00  after line 30 , it will work in this case too.

 Jn(x) is not accurately computed for large arguments, "JNX1" is better:

        b) Bessel Functions Jn(x) of integer order ( n >= 0 )
 

-The following program is based upon a program given by Keith Jarett in "HP-41 Extended Functions made easy"
-It uses the relations:   J0(x) + 2 J2(x) + 2 J4(x) + 2 J6(x) + ...... = 1  and   Jn-1(x) + Jn+1(x) = (2n/x) Jn(x)

Data Registers:  R00 thru R07 ( temp )  when the program stops, R01 = Jn(x) and R02 = x
Flags: /
Subroutines: /
 

01  LBL "JNX1"
02  STO 02
03  ABS
04  5
05  +
06  X<>Y
07  STO 00
08  X<Y?
09  X<>Y
10  INT
11  ST+ X
12  STO 06
13  1/X
14  STO 03
15  STO 04
16  STO 05
17   E40
18  STO 07
19  ISG 00
20  LBL 01
21  RCL 06
22  ST+ X
23  RCL 02
24  /
25  RCL 03
26  ST* Y
27  X<> 04
28  -
29  STO 03
30  RCL 06
31  2
32  MOD
33  *
34  ST+ 05
35  RCL 05
36  ABS
37  RCL 07                 lines 35 to 44 simply avoid an "OUT OF RANGE"
38  X>Y?
39  GTO 01
40  ST/ 01
41  ST/ 03
42  ST/ 04
43  ST/ 05
44  LBL 01
45  RCL 00
46  RCL 06
47  X#Y?
48  GTO 01
49  RCL 03
50  STO 01
51  LBL 01
52  DSE 06
53  GTO 01
54  RCL 05
55  ST+ X
56  RCL 03
57  -
58  ST/ 01
59  ST/ 03
60  ST/ 04
61  RCL 01
62  END

( 91 bytes / SIZE 008 )
 
 
      STACK        INPUTS     OUTPUTS
           Y             n           /
           X             x         Jn(x)

Example:    Calculate   J2(10)

     2    ENTER^
    10   XEQ "JNX1"   produces   0.254630314  (  in 25 seconds )

-Unlike "JNX" , "JNX1" yields accurate results for large x-values  ( for example, J0(1000) = 0.024786686 in 27mn26s )
 

       c) Bessel & modified Bessel Functions of the second kind  Yn(x) &  Kn(x)  ( non-integer order )
 

Formulae:    Yn(x) = ( Jn(x) cos(n(pi)) - J-n(x) ) / sin(n(pi))      ;      Kn(x) = (pi/2) ( I-n(x) - In(x) ) / sin(n(pi))      n # .... -3 ; -2 ; -1 ; 0 ; 1 ; 2 ; 3 ....
 

Data Registers:   R00 = x ; R01 = n ; R02 = temp
Flag: F01      If  F01 is clear,  "YNX" calculates  Yn(x)
                      If  F01 is set,     "YNX" ---------   Kn(x)
Subroutine: "JNX"
 
 

01  LBL "YNX"
02  DEG
03  STO 00
04  X<>Y
05  STO 01
06  CHS
07  X<>Y
08  XEQ "JNX"
09  STO 02
10  RCL 01
11  RCL 00
12  XEQ "JNX"
13  PI
14  R-D
15  RCL 01
16  *
17  STO Z
18  FS? 01
19  CLX
20  COS
21  *
22  RCL 02
23  -
24  X<>Y
25  SIN
26  /
27  PI
28  2
29  /
30  CHS
31  FC? 01
32  ST/ X
33  *
34  END

( 54 bytes / SIZE 003 )
 
 
      STACK        INPUTS      OUTPUTS
           Y             n             /
           X             x           f(x)

  with  f(x) =  Yn(x)  if  CF 01  ;   f(x) =  Kn(x)  if SF 01

Example:  Evaluate  Y1.4(3)  and   K1.4(3)

    CF 01   1.4  ENTER^   3   XEQ "YNX"   >>>>   Y1.4(3) = 0.137821836    ( in 31 seconds )
    SF 01    1.4  ENTER^   3          R/S          >>>>   K1.4(3) = 0.046088036      --------------
 

       d) Bessel & modified Bessel Functions of the second kind  Yn(x) &  Kn(x)  ( integer order )
 

Formulae:       with  psi(x) = the digamma function, we have ( positive n ):

     Yn(x) = -(1/pi) (x/2)-n SUMk=0,1,....,n-1  (n-k-1)!/(k!) (x2/4)k + (2/pi) ln(x/2) Jn(x) - (1/pi) (x/2)n SUMk=0,1,.....  ( psi(k+1) + psi(n+k+1) ) (-x2/4)k / (k!(n+k)!)

Kn(x) =  (1/2) (x/2)-n SUMk=0,1,....,n-1  (n-k-1)!/(k!) (-x2/4)k - (-1)n ln(x/2) In(x) + (1/2) (-1)n (x/2)n SUMk=0,1,.....  ( psi(k+1) + psi(n+k+1) ) (x2/4)k / (k!(n+k)!)
 
 

Data Registers:   R00 = x/2 ; R01 = n ; R02 to R05: temp
Flag: F01      If  F01 is clear,  "YNX1" calculates  Yn(x)
                      If  F01 is set,     "YNX1" ---------   Kn(x)
Subroutine: "JNX"
 

  01  LBL "YNX1"
  02  STO 00
  03  2
  04  ST/ 00
  05  X<> Z
  06  STO 01
  07  X<>Y
  08  XEQ "JNX"
  09  RCL 00
  10  X^2
  11  LN
  12  *
  13  1
  14  FS? 01
  15  CHS
  16  RCL 01
  17  Y^X
  18  *
  19  STO 02
  20  RCL 01
  21  STO 03
  22  CLX
  23  LBL 01
  24  RCL 01
  25  RCL 03
  26  -
  27  FACT
  28  RCL 03
  29  1
  30  -
  31  X<0?
  32  CLST
  33  FACT
  34  ST/ Y
  35  CLX
  36  RCL 00
  37  ST* X
  38  FS? 01
  39  CHS
  40  LASTX
  41  Y^X
  42  *
  43  +
  44  DSE 03
  45  GTO 01
  46  RCL 00
  47  RCL 01
  48  Y^X
  49  /
  50  ST- 02
  51  1.15443133
  52  STO 05
  53  RCL 01
  54  LBL 02
  55  X#0?
  56  1/X
  57  ST- 05
  58  X<> L
  59  DSE X
  60  GTO 02
  61  1
  62  STO 03
  63  STO 04
  64  RCL 05
  65  LBL 03
  66  RCL 04
  67  RCL 00
  68  X^2
  69  FC? 01
  70  CHS
  71  *
  72  RCL 03
  73  ST/ Y
  74  RCL 01
  75  +
  76  /
  77  STO 04
  78  RCL 05
  79  LASTX
  80  1/X
  81  -
  82  RCL 03
  83  1/X
  84  -
  85  STO 05
  86  *
  87  +
  88  ISG 03
  89  CLX
  90  X#Y?
  91  GTO 03
  92  RCL 01
  93  FACT
  94  /
  95  RCL 00
  96  FS? 01
  97  CHS
  98  RCL 01
  99  Y^X
100  *
101  RCL 02
102  +
103  FC? 01
104  PI
105  FS? 01
106  -2
107  /
108  END

( 151 bytes / SIZE 006 )
 
 
 
      STACK        INPUTS      OUTPUTS
           Y             n             /
           X             x           f(x)

  with  f(x) =  Yn(x)  if  CF 01  ;   f(x) =  Kn(x)  if SF 01   ;   n = 0 , 1 , 2 , ........ , 69

Example:  Evaluate  Y2(3)  and   K2(3)

    CF 01    2  ENTER^   3   XEQ "YNX1"   >>>>   Y2(3) = -0.160400393    ( in 23 seconds )
    SF 01    2  ENTER^   3          R/S            >>>>    K2(3) =  0.061510458      --------------

Note:    -If n is a positive integer, Y-n(x) = (-1)n Yn(x)  and  K-n(x) = Kn(x)
 

       e) An Asymptotic Expansion for  Jn(x) and Yn(x)

Formulae:     Jn(x) = (2/(pi*x))1/2 ( P.cos(x-(2n+1)pi/4) - Q.sin(x-(2n+1)pi/4) )  and   Yn(x) = (2/(pi*x))1/2 ( P.sin(x-(2n+1)pi/4) + Q.cos(x-(2n+1)pi/4) )
 

  where   P ~  1 - (4n2-1)(4n2-9)/(2!(8x)2) + (4n2-1)(4n2-9)(4n2-25)(4n2-49)/(4!(8x)4) - ......

     and   Q ~  (4n2-1)/(8x) - (4n2-1)(4n2-9)(4n2-25)/(3!(8x)3) + ......
 

Data Registers:  R00 = x ; R01 = n ; R02 to R06: temp
Flags: /
Subroutines: /
 

01  LBL "AEJY"
02  DEG
03  STO 00
04  X<>Y
05  STO 01
06  ST+ X
07  X^2
08  STO 04
09  1
10  STO 02
11  -
12  RCL 00
13  8
14  *
15  STO 05
16  ST* 05
17  /
18  STO 03
19  XEQ 01
20  STO 06
21  CLX
22  STO 02
23  SIGN
24  STO 03
25  XEQ 01
26  GTO 02
27  LBL 01
28  RCL 03
29  RCL 02
30  2
31  +
32  STO 02
33  ST+ X
34  3
35  -
36  X^2
37  RCL 04
38  -
39  *
40  RCL 04
41  RCL 02
42  ST+ X
43  DSE X
44  X^2
45  -
46  *
47  RCL 02
48  ST/ Y
49  DSE X
50  /
51  RCL 05
52  /
53  STO 03
54  +
55  X#Y?
56  GTO 01
57  RTN
58  LBL 02
59  RCL 00
60  RCL 01
61  ST+ X
62  1
63  +
64  4
65  /
66  PI
67  *
68  -
69  R-D
70  STO 04
71  X<>Y
72  P-R
73  RCL 04
74  RCL 06
75  P-R
76  ST+ T
77  RDN
78  -
79  2
80  PI
81  /
82  RCL 00
83  /
84  SQRT
85  ST* Z
86  *
87  END

( 112 bytes / SIZE 007 )
 
 
 
      STACK        INPUTS      OUTPUTS
           Y             n         Yn(x)
           X             x          Jn(x)

Example:       4   ENTER^  100   XEQ "AEJY"  >>>>   J4(100)  =  0.026105809
                                                             X<>Y                Y4(100)  = -0.075430120    ( in 10 seconds )
 

-The program will not work if  x is too "small" or if n is too "large"
 

       f) An Asymptotic Expansion for  Kn(x)

Formula:     Kn(x) ~  (pi/(2x))1/2 e-x  ( 1 + (4n2-1)/(8x) + (4n2-1)(4n2-9)/(2!(8x)2) + (4n2-1)(4n2-9)(4n2-25)/(3!(8x)3) + ......  )
 

Data Registers:  R00 = x ; R01 = 4n2 ; R02 , R03: temp
Flags: /
Subroutines: /
 

01  LBL "AEK"
02  STO 00
03  X<>Y
04  ST+ X
05  X^2
06  STO 01
07  SIGN
08  STO 02
09  STO 03
10  LBL 01
11  RCL 03
12  RCL 01
13  RCL 02
14  ST+ X
15  DSE X
16  X^2
17  -
18  *
19  RCL 00
20  RCL 02
21  *
22  8
23  *
24  /
25  STO 03
26  +
27  ISG 02
28  CLX
29  X#Y?
30  GTO 01
31  PI
32  RCL 00
33  ST+ X
34  /
35  SQRT
36  *
37  RCL 00
38  E^X
39  /
40  END

( 54 bytes / SIZE 004 )
 
 
 
      STACK        INPUTS      OUTPUTS
           Y             n             /
           X             x          Kn(x)

Example:       2   ENTER^  10   XEQ "AEK"  >>>>   K2(10)  = 0.000,021,509,817,05    ( in 15 seconds )
                      1.4  ENTER^  19       R/S          >>>>   K1.4(19) = 1.683198846 10-9          ( in 6 seconds )

-The program will not work if  x is too "small" or if n is too "large".
-For instance, n = 2 , x = 7 dont work.
-However, if we stop the iterations before the terms begin to increase, we find  K2(7) = 0.0005545622
-Divergence is more typical than convergence for asymptotic series...

Remark:     In(x) is accurately computed by "JNX" with flag F01 set, but if you need an asymptotic expansion for this function:

  replace line 39 with *
              line 36 with /
              line 34 with *  and add CHS after line 24   ( example:   I1.4(19) = 15,597,340  ( in 8 seconds ) )
 
 

       g) Integral of the Bessel Functions:   Jn(x)
 

Formula:        §0x Jn(t).dt = 2 ( Jn+1(x) + Jn+3(x) + Jn+5(x) + ........ )
 

Data Registers:   R00 = x/2 ; R01 = n ; R02 , R03: temp
Flag: F01
Subroutine: "JNX"
 

01  LBL "ITJ"
02  CF 01
03  STO 00
04  X<>Y
05  STO 01
06  1
07  +
08  STO 03
09  CLX
10  STO 02
11  LBL 01
12  RCL 03
13  RCL 00
14  XEQ "JNX"
15  2
16  ST+ 03
17  CLX
18  RCL 02
19  +
20  STO 02
21  LASTX
22  X#Y?
23  GTO 01
24  ST+ X
25  END

( 41 bytes / SIZE 004 )
 
 
 
      STACK        INPUTS      OUTPUTS
           Y             n             /
           X             x        §0x Jn(t).dt 

Example:     Calculate  A = §03  J1.4(x).dx

    1.4   ENTER^
     3     XEQ "ITJ"  >>>>  A = 1.049262784  ( in 1mn26s )
 

3°) Kummer's Function:  M(a;b;x)
 

Formula:      M(a;b;x) = 1 +  (a)1/(b)1. x1/1! + ............. +  (a)n/(b)n . xn/n! + ..........  where (a)n = a(a+1)(a+2) ...... (a+n-1)
 
 

        a) Real Variable
 

Data Registers:             R00:  x
                                      R01 = a
                                      R02 = b              registers R01 R02 are to be initialized before executing "KUM"
Flags: /
Subroutines:  /
 

01  LBL "KUM"
02  STO 00
03  CLST
04  SIGN
05  ENTER^
06  STO T
07  LBL 01
08  X<> T
09  RCL 01
10  R^
11  ST+ Y
12  RDN
13  *
14  RCL 02
15  R^
16  ST+ Y
17  ISG X
18  CLX
19  ST* Y
20  RDN
21  /
22  RCL 00
23  *
24  STO T
25  X<>Y
26  ST+ Y
27  X#Y?
28  GTO 01
29  END

( 46 bytes / SIZE 002 )
 
 
      STACK        INPUTS     OUTPUTS
           X             x       M(a;b;x)
           L             /            x

Example:     Compute  M(2;3;-Pi)

  2  STO 01
  3  STO 02
  PI   CHS   XEQ "KUM"  yields  0.166374562   ( in 13 seconds )

Notes:

 a)   2x (Pi)-1/2 M(1/2;3/2;-x2) = erf(x) = error function
 b)        (x/2)n  M(n+1/2;2n+1;2x) = Gamma(1+n) ex In(x)        where     In  = a modified Bessel function
 c)       (xa/a)   M(a;a+1;-x)  =  incgam(a;x)  = §0x e-t ta-1 dt                ( incgam = incomplete gamma function )

and many other functions are related to Kummer's functions.
 

        b) Complex Variable

-The parameters  a & b are still real, but the variable  z = x + i.y  is complex

Data Registers:        R00 and R03 thru R08: temp
                                      R01 = a
                                      R02 = b              registers R01 R02 are to be initialized before executing "KUMZ"
Flags: /
Subroutines:  /
 

01  LBL "KUMZ"
02  R-P
03  STO 00
04  X<>Y
05  STO 03
06  CLX
07  STO 05
08  STO 06
09  STO 08
10  SIGN
11  STO 04
12  STO 07
13  LBL 01
14  RCL 03
15  RCL 08
16  +
17  STO 08
18  RCL 01
19  RCL 02
20  RCL 06
21  ST+ Z
22  +
23  ISG 06
24  CLX
25  RCL 06
26  *
27  /
28  RCL 00
29  *
30  RCL 07
31  *
32  STO 07
33  P-R
34  RCL 04
35  +
36  STO04
37  LASTX
38  -
39  X^2
40  X<>Y
41  RCL 05
42  +
43  STO 05
44  LASTX
45  -
46  X^2
47  +
48  X#0?
49  GTO 01
50  RCL 05
51  RCL 04
52  END

( 64 bytes / SIZE 009 )
 
 
      STACK        INPUTS      OUTPUTS
           Y             y             y'
           X             x             x'

   with  Kum ( a ; b ; x+i.y ) = x' + i.y'

Example:     If  a = 4 ;  b = 3     4  STO 01   3  STO 02

    2  ENTER^
    1  XEQ "KUMZ  >>>>   -3.156090293   X<>Y   2.541499313

 Whence    Kum ( 4 ; 3 ; 1 + 2.i ) =  -3.156090293 + i. 2.541499313
 

 4°) Parabolic Cylinder Functions:  U(a;x)  V(a;x)  W(a;x)

        a) U(a;x) & V(a;x)

-The standard solutions of the differential equation   d2y/dx2 - ( x2/4 + a ).y = 0  may be expressed with the Kummer's Functions

     U(a;x) = pi1/22-1/4-a/2e-x^2/4 (Gam(3/4+a/2))-1 M(1/4+a/2;1/2;x2/2) - pi1/221/4-a/2 x.e-x^2/4 (Gam(1/4+a/2))-1 M(3/4+a/2;3/2;x2/2)

     V(a;x) =  Gam(a+1/2) [ sin(a.pi) U(a;x) + U(a;-x) ].(pi)-1

-However, if  a = -1/2 ; -3/2 ; ....  the Gamma function will be infinite. Therefore, other formulas are used if a < 0 , namely:

     U(a;x) = Y1 cos(pi(1/4+a/2)) - Y2 sin(pi(1/4+a/2))     and     V(a;x) = (Gam(1/2-a))-1 [ Y1 sin(pi(1/4+a/2)) + Y2 cos(pi(1/4+a/2)) ]

  where  Y1 = pi-1/22-1/4-a/2e-x^2/4 (Gam(1/4-a/2)) M(1/4+a/2;1/2;x2/2)  and   Y2 = pi-1/221/4-a/2 x.e-x^2/4 (Gam(3/4-a/2)) M(3/4+a/2;3/2;x2/2)
 
 

Data Registers:  R00 thru R06  ( R03 = x ; R04 = a )
Flags: /
Subroutines:  "GAM"  ( Gamma Function )
                        "KUM"  ( Kummer's Functions )
 

001  LBL "PCF1"
002  STO 03
003  X^2
004  X<>Y
005  STO 04
006  .5
007  STO 02
008  ST* Y
009  ST* Z
010  X^2
011  LASTX
012  +
013  +
014  STO 01
015  X<>Y
016  ISG 02
017  XEQ "KUM"
018  STO 05
019  .5
020  ST- 01
021  STO 02
022  RCL 00
023  XEQ "KUM"
024  RCL 00
025  RCL 02
026  *
027  CHS
028  E^X
029  ST* 05
030  *
031  2
032  RCL 01
033  Y^X
034  /
035  STO 06
036  RCL 03
037  LASTX
038  RCL 02
039  SQRT
040  *
041  /
042  ST* 05
043  RCL 02
044  RCL 01
045  RCL 04
046  SIGN
047  *
048  +
049  STO 00
050  XEQ "GAM"
051  PI
052  SQRT
053  /
054  RCL 04
055  SIGN
056  Y^X
057  ST/ 06
058  RCL 00
059  RCL 02
060  LASTX
061  *
062  -
063  XEQ "GAM"
064  PI
065  SQRT
066  /
067  RCL 04
068  SIGN
069  Y^X
070  ST/ 05
071  RCL 02
072  RCL 04
073  ABS
074  +
075  XEQ "GAM"
076  STO 02
077  1
078  CHS
079  ACOS
080  STO 00
081  ST* 01
082  RCL 04
083  X<0?
084  GTO 01
085  RCL 06
086  RCL 05
087  -
088  STO Y
089  RCL 00
090  RCL 04
091  *
092  SIN
093  *
094  RCL 05
095  RCL 06
096  +
097  +
098  RCL 02
099  *
100  PI
101  /
102  GTO 02
103  LBL 01
104  RCL 01
105  RCL 05
106  P-R
107  RCL 01
108  RCL 06
109  P-R
110  R^
111  -
112  RDN
113  +
114  RCL 02
115  /
116  LBL 02
117  R^
118  END
 

( 161 bytes / SIZE 007 )
 
 
      STACK        INPUTS      OUTPUTS
           Y             a         V(a;x)
           X             x         U(a;x)

Example1:     Calculate  U(0.4;1.9) & V(0.4;1.9)

    0.4   ENTER^
    1.9   XEQ "PCF1"  >>>>  U(0.4;1.9) = 0.194020564     X<>Y      V(0.4;1.9) =  1.882850363    ( in 42 seconds )

Example2:     Calculate  U(-0.4;1.9) & V(-0.4;1.9)

    -0.4   ENTER^
     1.9   XEQ "PCF1"  >>>>  U(-0.4;1.9) = 0.376027811     X<>Y     V(-0.4;1.9) =  1.376169516  ( in 41 seconds )
 
 

        b) W(a;x)

-The differential equation   d2y/dx2 + ( x2/4 - a ).y = 0  is solved by a series expansion   y = a0 + a1.x + a2.x2 + ...... + ak.xk + ......
-The standard solution  W(a;x) is defined by

       a0 = 2-3/4 | Gam(1/4 + i.a/2) / Gam(3/4 + i.a/2) |1/2   and   a1 = -2-1/4 | Gam(3/4 + i.a/2) / Gam(1/4 + i.a/2) |1/2

-The relation  | Gam(1/4 + i.a/2) |.| Gam(3/4 + i.a/2) | =  pi.(2/cosh(a.pi))1/2  is also used to call "GAMZ" only once.

Data Registers:  R00 thru R09 ( temp )  When the program stops, R01 = W(a;x)
Flags: /
Subroutines:  "GAMZ"  ( Gamma function in the complex plane )
 

01  LBL "PCF2"
02  STO 06
03  X<>Y
04  STO 07
05  .5
06  ST* Y
07  X^2
08  XEQ "GAMZ"
09  LASTX
10  PI
11  RCL 07
12  *
13  E^X
14  ENTER^
15  1/X
16  +
17  .5
18  STO 04
19  *
20  SQRT
21  PI
22  /
23  SQRT
24  *
25  ST* 04
26  1/X
27  CHS
28  STO 05
29  RCL 06
30  STO 08
31  *
32  RCL 04
33  +
34  STO 01
35  CLX
36  STO 02
37  STO 03
38  SIGN
39  STO 00
40  LBL 01
41  3
42  STO 09
43  LBL 02
44  RCL 04
45  RCL 07
46  *
47  RCL 02
48  4
49  /
50  -
51  RCL 00
52  RCL 00
53  1
54  +
55  STO 00
56  *
57  /
58  X<> 05
59  X<> 04
60  X<> 03
61  STO 02
62  RCL 06
63  RCL 08
64  *
65  STO 08
66  RCL 05
67  *
68  RCL 01
69  +
70  STO 01
71  LASTX
72  X#Y?
73  GTO 01
74  DSE 09
75  GTO 02
76  END

( 100 bytes / SIZE 010 )
 
 
 
      STACK        INPUTS     OUTPUTS
           Y             a        W(a;x)
           X             x        W(a;x)

Example:     Calculate  W(0.4;1.9)

    0.4   ENTER^
    1.9   XEQ "PCF2"  >>>>   0.219336459   ( in 52 seconds )
 

       c) Asymptotic Expansions for  U(a;x) & V(a;x)   ( x large , a moderate )
 

Formulas:    U(a;x)  ~   e-x^2/4  x-a-1/2 [ 1 - (a+1/2)(a+3/2)/(1.(2x2)) + (a+1/2)(a+3/2)(a+5/2)(a+7/2)/(2.(2x2)2) -  ....... ]

                   V(a;x)  ~   (2/pi)1/2  ex^2/4  xa-1/2 [ 1 + (a-1/2)(a-3/2)/(1.(2x2)) + (a-1/2)(a-3/2)(a-5/2)(a-7/2)/(2.(2x2)2) + ....... ]
 

Data Registers:  R00 = x ; R01 = a ; R02 , R03: temp
Flags:  F01    if F01 is clear, AEUV computes  U(a;x)
                       if F01 is set,    -----------------  V(a;x)
Subroutines: /
 

01  LBL "AEUV"
02  STO 00
03  X<>Y
04  STO 01
05  1
06  STO 02
07  STO 03
08  LBL 01
09  RCL 03
10  ST+ X
11  ENTER^
12  DSE X
13  .5
14  ST- Z
15  -
16  RCL 01
17  FS? 01
18  CHS
19  ST+ Z
20  +
21  *
22  RCL 02
23  FC? 01
24  CHS
25  *
26  RCL 03
27  /
28  RCL 00
29  X^2
30  ST+ X
31  /
32  ISG 03
33  CLX
34  STO 02
35  +
36  X#Y?
37  GTO 01
38  RCL 00
39  X^2
40  FC? 01
41  CHS
42  4
43  /
44  E^X
45  *
46  RCL 00
47  RCL 01
48  FC? 01
49  CHS
50  .5
51  -
52  Y^X
53  *
54  2
55  PI
56  /
57  SQRT
58  FC? 01
59  SIGN
60  *
61  END

( 84 bytes / SIZE 004 )
 
 
      STACK        INPUTS      OUTPUTS
           Y             a             /
           X             x           f(x)

where  f(x) = U(a;x)  if  flag F01 is clear
           f(x) = V(a;x)  if  flag F01 is set.
 

Example:    CF 01   2  ENTER^   10   XEQ "AEUV"  >>>>   U(2;10) = 4.210624069 10-14   ( in 14 seconds )
                    SF 01   2  ENTER^   10         R/S            >>>>   V(2;10) = 1.823604920 1012      ( in 7 seconds )
 

Note:     V(2;10) is accurately computed by "PCF1" ( it yields  1.823604925 1012 in 2mn26s )
           but "PCF1" gives  U(2;10) = -2000 ( a completely wrong result! )
     This meaningless value results from the subtraction of 2 large numbers, and the leading digits are cancelled.
-If x is too small, the series will diverge too soon.
-However,  U(-5;5) = 1.879976816  is correctly calculated ( the negative a-value helps convergence).
 

       d) An Asymptotic Expansion for  W(a;x)    ( x large , a moderate )
 

Formulae:      W(a;x)  ~  (2k/x)1/2  ( s1(a;x) cos A - s2(a;x) sin A )      and     W(a;-x)  ~  (2/(kx))1/2 ( s1(a;x) sin A + s2(a;x) cos A )      ( x > 0 )

          where        A = x2/2 - a.Ln(x) + pi/4 + B/2   with  B = arg ( Gam(1/2+i.a) )
                            k = 1/(epi.a+(1+e2pi.a)1/2)
     s1(a;x)+i.s2(a;x)  =  Sumn=0,1,2,......  (-i)n Gam(2n+1/2+i.a)/Gam(1/2+i.a)  1/(n!(2x2)n)
 

Data Registers:  R00 thru R08: temp  when the program stops, R02 = a ; R08 = x
Flags: /
Subroutine:  "GAMZ"  ( gamma function in the complex plane )
 
 

  01  LBL "AEW"
  02  STO 08
  03  X<>Y
  04  STO 02
  05  CLX
  06  STO 01
  07  STO 04
  08  STO 07
  09  SIGN
  10  STO 03
  11  STO 05
  12  STO 06
  13  LBL 01
  14  ISG 01
  15  CLX
  16  RCL 02
  17  RCL 01
  18  ST+ X
  19  .5
  20  -
  21  R-P
  22  1
  23  LASTX
  24  -
  25  RCL 02
  26  R-P
  27  ST* Z
  28  X<> T
  29  +
  30  RCL 04
  31  +
  32  STO 04
  33  X<>Y
  34  RCL 03
  35  *
  36  STO 03
  37  RCL 05
  38  RCL 08
  39  X^2
  40  ST+ X
  41  *
  42  RCL 01
  43  *
  44  STO 05
  45  /
  46  P-R
  47  RCL 06
  48  +
  49  STO 06
  50  LASTX
  51  -
  52  ABS
  53  X<>Y
  54  RCL 07
  55  +
  56  STO 07
  57  LASTX
  58  -
  59  ABS
  60  +
  61  X#0?
  62  GTO 01
  63  RCL 02
  64  .5
  65  XEQ "GAMZ"
  66  R-P
  67  CLX
  68  2
  69  /
  70  RCL 08
  71  X^2
  72  PI
  73  +
  74  4
  75  /
  76  RCL 02
  77  RCL 08
  78  LN
  79  *
  80  -
  81  R-D
  82  +
  83  RCL 08
  84  SIGN
  85  X<0?
  86  CLX
  87  ASIN
  88  +
  89  1
  90  P-R
  91  RCL 07
  92  *
  93  X<>Y
  94  RCL 06
  95  *
  96  +
  97  RCL 02
  98  PI
  99  *
100  E^X
101  ENTER^
102  X^2
103  1
104  +
105  SQRT
106  +
107  RCL 08
108  SIGN
109  CHS
110  Y^X
111  ST+ X
112  RCL 08
113  /
114  SQRT
115  *
116  END

( 138 bytes / SIZE 009 )
 
 
      STACK        INPUTS      OUTPUTS
           Y             a             /
           X             x         W(a;x)

Example:     0.5  ENTER^   10   XEQ "AEW"  >>>>   0.092208658    ( in 59 seconds )
                    -0.5  ENTER^   10          R/S         >>>>  -0.228640282     --------------
 

-The series diverge too soon if x is too small, but if you add   RND  after line 60

   FIX 5    1  ENTER^   5   XEQ "AEW"  yields   W(1;5) = 0.022808   and similarly   W(-1;5) = -0.570255   ( in 57 seconds )
 

        e) A Recurrence Relation for U(a;x)
 

-If a and x are large, the recurrence relation    (a+1/2) U(a+1;x) + x.U(a;x) = U(a-1;x)   may be used straightforward from 2 values.
-However, this can lead to low accuracy by cancellation of leading digits.
-Better accuracy is attainable if we use this relation in the reverse direction, starting with two arbitrary values ( 1 and 0 ) for a+21 and a+20 ( line 05 )
-Then, "AEUV" is called as a subroutine and finally, a normalization factor gives U(a;x)
 

Data Registers:  R00 = x  ; R04 = a ; R02 , R03 , R07: temp
Flag: CF01
Subroutine:  "AEUV"
 

01  LBL "UAX"
02  STO 00
03  X<>Y
04  STO 04
05  20
06  +
07  STO 01
08  CLX
09  STO 02
10  SIGN
11  STO 03
12  LBL 01
13  RCL 02
14  RCL 01
15  .5
16  +
17  *
18  RCL 03
19  STO 02
20  RCL 00
21  *
22  +
23  STO 03
24  RCL 01
25  1
26  -
27  STO 01
28  RCL 04
29  X#Y?
30  GTO 02
31  RCL 03
32  STO 07
33  LBL 02
34  RCL 03
35  ABS
36   E10
37  X>Y?
38  GTO 01
39  RCL 00
40  ST+ X
41  RCL 01
42  +
43  X>0?
44  GTO 01                 the recurrence relation is applied until  2x+a' <= 0 but other choices may be better ( for instance | a | < 1 )
45  RCL 04                 especially if line 54 is replaced by  XEQ "PCF1"
46  RCL 01
47  X>Y?
48  GTO 01
49  RCL 03
50  ST/ 07
51  RCL 01
52  RCL 00
53  CF 01
54  XEQ "AEUV"
55  RCL 07
56  *
57  END

( 81 bytes / SIZE 008 )
 
 
      STACK        INPUTS      OUTPUTS
           Y             a             /
           X             x         U(a;x)

Examples:    5   ENTER^   XEQ "UAX"   >>>>   U(5;5) =  1.552271290 10-7       ( execution time = 45 seconds )
                    12  ENTER^   7     R/S          >>>>   U(12;7) = 3.282492495 10-17     ( execution time = 55 seconds )
 

-In these examples, "PCF1" would yield meaningless results for U(a;x) because of cancellation of the leading digits.
-If  x is too "small" , the asymptotic expansion used in "AEUV" may diverge too soon.
-In this case, replace line 54 with  XEQ "PCF1"  ( that's why I've used register R07 instead of R05 )
 
 

  5°) Hypergeometric Functions:  F(a;b;c;x)
 

        a) Program#1
 

Formula:      F(a;b;c;x) = 1 +  (a)1(b)1/(c)1. x1/1! + ............. +  (a)n(b)n/(c)n . xn/n! + ..........  where (a)n = a(a+1)(a+2) ...... (a+n-1)    and     | x | < 1
 

Data Registers:             R00:  x
                                      R01 = a
                                      R02 = b              registers R01 R02 R03 are to be initialized before executing "HGF"
                                      R03 = c

Flags: /
Subroutines:  /

01  LBL "HGF"
02  STO 00
03  CLST
04  SIGN
05  ENTER^
06  STO T
07  LBL 01
08  X<> T
09  RCL 01
10  R^
11  ST+ Y
12  RDN
13  *
14  RCL 02
15  R^
16  ST+ Y
17  RDN
18  *
19  RCL 03
20  R^
21  ST+ Y
22  ISG X
23  CLX
24  ST* Y
25  RDN
26  /
27  RCL 00
28  *
29  STO T
30  X<>Y
31  ST+ Y
32  X#Y?
33  GTO 01
34  END

( 52 bytes / SIZE 004 )
 
 
      STACK        INPUTS      OUTPUTS
           X             x       F(a;b;c;x)
           L             /             x

Example:   Calculate  F(0.3 ; -0.7 ; 0.4 ; 0.49)

  0.3   STO 01    -0.7   STO 02    0.4    STO 03
  0.49   XEQ "HGF"  yields  0.720164356     ( in 18 seconds )
 

        b) Program#2

-This second program uses the following properties:

- If    c-a-b > 0  ,  F(a;b;c;1) = Gam(c).Gam(c-a-b)/(Gam(c-a).Gam(c-b))
- F(a;b;b;x) = F(b;a;b;x) = 1/(1-x)a
- If  c is a negative integer and neither a nor b is a negative integer such that  -a<-c , -b<-c
  then F is not defined but "HGF2" sets flag F00 and calculates:

         Limt tends to c ( F(a,b,t,x) )/Gam(t)  = (a)1-c(b)1-c/(1-c)!  F(a-c+1;b-c+1;2-c;x)

-This may be used in some of the Associated Legendre Function programs hereafter
 

Data Registers:             R00:  x
                                      R01 = a
                                      R02 = b              registers R01 R02 R03 are to be initialized before executing "HGF"
                                      R03 = c

Flags:  F00 & F05
Subroutine:  "GAM"  ( only if  x = 1 )
 

  01  LBL "HGF2"
  02  STO 00
  03  CF 00
  04  CF 05
  05  RCL 01
  06  X>0?
  07  GTO 00
  08  FRC
  09  X#0?
  10  GTO 00
  11  SF 05
  12  RCL 03
  13  LASTX
  14  -
  15  X<0?
  16  GTO 03
  17  LBL 00
  18  RCL 02
  19  X>0?
  20  GTO 01
  21  FRC
  22  X#0?
  23  GTO 01
  24  SF 05
  25  RCL 03
  26  LASTX
  27  -
  28  X<0?
  29  GTO 03
  30  LBL 01
  31  RCL 03
  32  X>0?
  33  GTO 02
  34  FRC
  35  X#0?
  36  GTO 02
  37  SF 00
  38  CF 05
  39  1
  40  RCL 03
  41  -
  42  ST+ 01
  43  ST+ 02
  44  1
  45  +
  46  STO 03
  47  LBL 02
  48  RCL 00
  49  1
  50  X=Y?
  51  FS? 05
  52  GTO 03
  53  RCL 03
  54  XEQ "GAM"
  55  STO N
  56  RCL 03
  57  RCL 01
  58  RCL 02
  59  +
  60  -
  61  XEQ "GAM"
  62  ST* N
  63  RCL 03
  64  RCL 01
  65  -
  66  XEQ "GAM"
  67  ST/ N
  68  RCL 03
  69  RCL 02
  70  -
  71  XEQ "GAM"
  72  ST/ N
  73  0
  74  X<> N
  75  GTO 05
  76  LBL 03
  77  RCL 01
  78  RCL 02
  79  RCL 03
  80  X=Y?
  81  GTO 04
  82  RDN
  83  X<>Y
  84  R^
  85  X#Y?
  86  GTO 06
  87  LBL 04
  88  X<> Z
  89  1
  90  RCL 00
  91  -
  92  X<>Y
  93  CHS
  94  Y^X
  95  LBL 05
  96  FS? 00
  97  GTO 08
  98  GTO 10
  99  LBL 06
100  CLST               ( lines 100 to 130 may be replaced by  RCL 00   XEQ "HGF" )
101  SIGN
102  ENTER^
103  STO T
104  LBL 07
105  X<> T
106  RCL 01
107  R^
108  ST+ Y
109  RDN
110  *
111  RCL 02
112  R^
113  ST+ Y
114  RDN
115  *
116  RCL 03
117  R^
118  ST+ Y
119  ISG X
120  CLX
121  ST* Y
122  RDN
123  /
124  RCL 00
125  *
126  STO T
127  X<>Y
128  ST+ Y
129  X#Y?
130  GTO 07
131  FC? 00
132  GTO 10
133  LBL 08
134  2
135  RCL 03
136  -
137  STO 03
138  1
139  -
140  ST+ 01
141  ST+ 02
142  CLX
143  RCL 00
144  LASTX
145  RCL 03
146  -
147  STO M
148  Y^X
149  *
150  LBL 09
151  RCL 01
152  RCL M
153  1
154  -
155  +
156  *
157  RCL 02
158  RCL M
159  ST/ Z
160  1
161  -
162  +
163  *
164  DSE M
165  GTO 09
166  LBL 10
167  CF 05
168  END

( 244 bytes / SIZE 004 )
 
 
      STACK        INPUTS      OUTPUTS
           X             x            f(x)

 where   f(x) = F(a;b;c;x)                                            if F00 is clear
             f(x) =  Limt tends to c ( F(a,b,t,x) )/Gam(t)      if F00 is set
 

Examples:    1  STO 01   2  STO 02   3  STO 03         0.1  XEQ "HGF2"  >>>>  F(1;2;3;0.1) = 1.072103131  ( in 7 seconds )

and similarly                                F(1;2;7;1)   = 1.5                    ( 2.5 s )
                                                  F(2;3;3;0.7) = 11.11111111   ( 1.5 s )
          Lim t tends to -7  ( F(2;2;t;0.1)/Gam(t) ) = 0.105229334   ( 19 s )          ( flag F00 is set )
          Lim t tends to -1  ( F(-4;-5;t;1)/Gam(t) ) = 420                  ( 5 s )            ( flag F00 is set )    .....  etc  .....
 
 
 

6°) Legendre Functions

        a) 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.

        b) Associated Legendre Functions Pnm(x) & Qnm(x) with 0 <= n <= m  (integer order & degree)
 

Formulae:     Pnm(x) = ( x2-1 )m/2 dmPn(x)/dxm  if  | x | > 1           Pnm(x) = (-1)m ( 1- x2 )m/2 dmPn(x)/dxm  if  | x | < 1
                      Qnm(x) = ( x2-1 )m/2 dmQn(x)/dxm  if  | x | > 1         Qnm(x) = (-1)m ( 1- x2 )m/2 dmQn(x)/dxm  if  | x | < 1

-if  m = 0   Pnm(x) = Pn(x) = Legendre polynomials
        and   Qnm(x) = Qn(x)  with  Q0(x) = 0.5 Ln | (1+x)/(1-x) |  ;  Q1(x) = x/2 Ln | (1+x)/(1-x) | - 1  and   n.Qn(x) = (2n-1).x.Qn-1(x) - (k-1).Qn-2(x)

-We have employed the recurrence relations:   Pnm+1(x) = | x2-1 |-1/2 [ (n-m).x. Pnm(x) - (n+m). Pn-1m(x) ]   ( the same relation holds for Qnm(x) )

Important note:   If | x | is significantly greater than 1 , Qnm(x) is obtained ( very ) inaccurately and "ALF2" is preferable ( see d) below )
                             ( the recurrende relation is unstable in this case )

Data Registers:      R00 to Rnn: temp
Flags: F05 & F02    if flag F02 is clear "ALF" calculates  Pnm(x)
                                 if flag F02 is set    "ALF" ----------  Qnm(x)
Subroutines: /
 

-This program uses 4 synthetic registers M N O P
-Don't interrupt "ALF": the content of register P could be altered
-These 4 registers may of course be replaced by any unused data registers ( Rpp with p > n )
 

001  LBL "ALF"
002  X<> Z
003  STO O
004  RDN
005  CF 05
006  X=0?
007  SF 05
008   E3
009  /
010  1
011  +
012  STO M
013  X<>Y
014  STO Y
015  LASTX
016  STO 00
017  FC? 02
018  GTO 01
019  ST+ Y
020  RCL Z
021  -
022  /
023  ABS
024  SQRT
025  LN
026  STO 00
027  LBL 01
028  STO N
029  RCL Z
030  *
031  STO P
032  -
033  RCL M
034  INT
035  1/X
036  ENTER^
037  SIGN
038  -
039  FS? 02
040  X#0?
041  GTO 02
042  SIGN
043  STO Y
044  CHS
045  LBL 02
046  *
047  RCL P
048  +
049  RCL N
050  X<>Y
051  STO IND M
052  ISG M
053  GTO 01
054  RCL Z
055  1
056  ST- M
057  X<> O
058  X=0?
059  GTO 04
060  1
061  -
062   E3
063  /
064  STO O
065  X<> M
066  INT
067  STO N
068  STO P
069  X<>Y
070  LBL 03
071  RCL IND N
072  RCL Y
073  *
074  RCL N
075  RCL M
076  INT
077  -
078  ST* Y
079  LASTX
080  ST+ X
081  +
082  DSE N
083  ""                       TEXT0   or another NOP instruction like  LBL 00  STO X  ...  etc  ...
084  RCL IND N
085  *
086  -
087  X<>Y
088  X^2
089  ENTER^
090  SIGN
091  ST+ N
092  -
093  ABS
094  SQRT
095  X=0?
096  SIGN
097  /
098  STO IND N
099  DSE N
100  ""                       TEXT0   or another NOP instruction like  LBL 00  STO X  ...  etc  ...
101  X<>Y
102  ISG O
103  GTO 03
104  RCL P
105  STO N
106  SIGN
107  RCL O
108  FRC
109  +
110  RCL M
111  INT
112  +
113  STO O
114  X<>Y
115  ISG M
116  GTO 03
117  X<>Y
118  LBL 04
119  RDN
120  SIGN
121  RDN
122  FS?C 05
123  RCL 00
124  CLA
125  END

( 186 bytes / SIZE nnn+1 )
 
 
 
      STACK        INPUTS      OUTPUTS
           Z             m            /
           Y             n             /
           X             x          F(x)
           L             /            x

with  F(x) =  Pnm(x)  if flag F02 is clear  and   F(x) =   Qnm(x)  if flag F02 is set.
 

Example:    Calculate   P74(0.6)  &   Q74(0.6)
 

a) SIZE 008  ( or greater )

b) CF 02    4      ENTER^
                  7      ENTER^
                0.6     XEQ "ALF"  yields   P74(0.6) = 715.3090560   ( in 17 seconds )

c) SF 02    4      ENTER^
                 7      ENTER^
               0.6       R/S          gives     Q74(0.6) = -1011.171802  ( in 17 seconds )

-Similarly:      P74(1.2) = 6327.212011   &   Q74(1.2) = 82.12107441

-For x = 1.2   Q74(x)  is still acceptable ( "ALF2" produces  82.12107808 )
  but for x = 3 , the errors are too great and "ALF2" is much better ( cf d) below )
 

        c) Associated Legendre Functions of the first kind ( fractional order & degree ) with  -1 < x  < 3
 

Formula:      Pnm(x) = |(x+1)/(x-1)|m/2  F(-n ; n+1 ; 1-m ; (1-x)/2 ) / Gamma(1-m)        (  x # 1 )
 

Data Registers:         R00 to R06: temp
Flags: /
Subroutines:  "GAM"  ( Gamma function )
                         "HGF"   ( Hypergeometric functions )
 

01  LBL "ALF1"
02  STO 04
03  CLX
04  SIGN
05  STO T
06  X<>Y
07  CHS
08  STO 01
09  -
10  STO 02
11  RDN
12  STO 05
13  -
14  STO 03
15  CLX
16  RCL 04
17  -
18  2
19  /
20  XEQ "HGF"
21  STO 06
22  RCL 03
23  XEQ "GAM"
24  ST/ 06
25  RCL 06
26  RCL 04
27  1
28  +
29  RCL 04
30  LASTX
31  -
32  /
33  ABS
34  RCL 05
35  2
36  /
37  Y^X
38  *
39  END

( 58 bytes / SIZE 007 )
 
 
 
      STACK        INPUTS      OUTPUTS
           Z             m             /
           Y             n             /
           X             x         Pnm(x)

Example:    Calculate   P1.30.4(1.9)

    0.4  ENTER^
    1.3  ENTER^
    1.9  XEQ "ALF1"  gives   2.980130385  ( in 26 seconds )

-This program will not work if m = 1,2,3,4, ..... unless you replace line 20 by  XEQ "HGF2"  and line 23 by  FC? 00  XEQ "GAM"  FC?C 00

-With this modification, you'll find, for example:   P73(1.7) = 102985.1588   ( in 9 seconds )
 

        d) Associated Legendre Functions of the second kind ( fractional order & degree ) with  | x | > 1
 

Formula:      Qnm(x) = ei.mpi 2-n-1 pi1/2 x-n-m-1 (x2-1)m/2 F(1+n/2+m/2 ; 1/2+n/2+m/2 ; n+3/2 ; 1/x2 )  Gamma(1+n+m) / Gamma(n+3/2)

    ( the result is a complex number )

-If  x < -1 and  -n-m-1 is not an integer, there will be a "DATA ERROR" message.
 

Data Registers:         R00 to R07: temp
Flags: /
Subroutines:  "GAM"  ( Gamma function )
                        "HGF"   ( Hypergeometric functions )
 

01  LBL "ALF2"
02  STO 04
03  CLX
04  SIGN
05  +
06  STO 03
07  STO 05
08  X<>Y
09  STO 06
10  +
11  .5
12  ST+ 03
13  *
14  STO 02
15  LASTX
16  +
17  STO 01
18  RCL 04
19  X^2
20  1/X
21  XEQ "HGF"
22  STO 07
23  RCL 05
24  RCL 06
25  +
26  XEQ "GAM"
27  ST* 07
28  RCL 03
29  XEQ "GAM"
30  ST/ 07
31  RCL 07
32  RCL 04
33  X^2
34  1
35  -
36  RCL 06
37  2
38  /
39  Y^X
40  *
41  RCL 04
42  RCL 05
43  RCL 06
44  +
45  Y^X
46  /
47  2
48  RCL 05
49  Y^X
50  /
51  PI
52  SQRT
53  *
54  1
55  CHS
56  ACOS
57  RCL 06
58  *
59  X<>Y
60  P-R
61  END

( 86 bytes / SIZE 008 )
 
 
 
      STACK        INPUTS      OUTPUTS
           Z             m             /
           Y             n            B
           X             x            A

                     with   Qnm(x) =  A + i.B
 

Examples:    Compute    Q1.20.7(1.9)

    0.7   ENTER^
    1.2   ENTER^
    1.9   XEQ "ALF2"  yields    -0.081513570
                                  X<>Y    0.112193804     thus,     Q1.20.7(1.9) = -0.081513570 + 0.112193804 . i      ( in 29 seconds )

-Similarly, we find  Q74(3) = 0.004063232  ( in 18 seconds )  ( "ALF" would give a wrong result! )

-This program will not work if  n+3/2 = 0,-1,-2,-3,-4, ..... unless you replace line 21 by  XEQ "HGF2"  and line 29 by  FC? 00  XEQ "GAM"  FC?C 00

-With this modification, you'll find, for example:   Q -1.52(7) = 0.114719166   ( in 16 seconds )

-If you need  | Qnm(x) | only, delete lines 54 to 60 and for instance:   |  Q1.20.7(1.9) | = 0.138679169   ( this real result is in L-register otherwise )
 

        e) Associated Legendre Functions of the first kind ( fractional order & degree )  with   x > 1
 

Formula:      Pnm(x) = 2-n-1pi-1/2Gam(-1/2-n) x-n+m-1/ ((x2-1)m/2Gam(-n-m)) F(1/2+n/2-m/2;1+n/2-m/2;n+3/2;1/x2)
                                    +2npi-1/2Gam(1/2+n) xn+m/ ((x2-1)m/2Gam(1+n-m)) F(-n/2-m/2;1/2-n/2-m/2;1/2-n;1/x2)

( The factor pi-1/2 in the second term is omitted by mistake in the "Handbook of Mathematical Functions" page 332 )
 

Data Registers:         R00 to R08: temp
Flags: /
Subroutines:  "GAM"  ( Gamma function )
                         "HGF"   ( Hypergeometric functions )

01  LBL "ALF3"
02  STO 04
03  RDN
04  STO 05
05  X<>Y
06  STO 06
07  -
08  2
09  ST/ Y
10  1/X
11  +
12  STO 01
13  LASTX
14  +
15  STO 02
16  1.5
17  RCL 05
18  +
19  STO 03
20  RCL 04
21  X^2
22  1/X
23  XEQ "HGF"
24  RCL 04
25  RCL 01
26  ST+ X
27  Y^X
28  ST+ X
29  /
30  STO 07
31  1
32  RCL 03
33  -
34  XEQ "GAM"
35  ST* 07
36  RCL 05
37  RCL 06
38  +
39  CHS
40  STO 01
41  XEQ "GAM"
42  ST/ 07
43  2
44  1/X
45  ST* 01
46  RCL 01
47  +
48  STO 02
49  .5
50  RCL 05
51  -
52  STO 03
53  RCL 00
54  XEQ "HGF"
55  2
56  RCL 05
57  Y^X
58  ST/ 07
59  *
60  RCL 04
61  RCL 05
62  RCL 06
63  +
64  Y^X
65  *
66  STO 08
67  RCL 05
68  RCL 06
69  -
70  1
71  +
72  XEQ "GAM"
73  ST/ 08
74  1
75  RCL 03
76  -
77  XEQ "GAM"
78  RCL 08
79  *
80  RCL 07
81  +
82  PI
83  SQRT
84  /
85  RCL 04
86  X^2
87  1
88  -
89  RCL 06
90  2
91  /
92  Y^X
93  /
94  END

( 138 bytes / SIZE 009 )
 
 
 
      STACK        INPUTS      OUTPUTS
           Z             m             /
           Y             n             /
           X             x         Pnm(x)

Example:    Calculate   P1.7-0.6(4.8)

  -0.6  ENTER^
   1.7  ENTER^
   4.8  XEQ "ALF3"  >>>>>  10.67810283    ( in 38 seconds )

-Many other relations between the Legendre Functions and the hypergeometric Functions may be found  in the "Handbook of Mathematical Functions".
 

        f) Associated Legendre Functions Pnm(x) & Qnm(x)  (fractional order & degree)  | x | < 1

-Actually, these functions are solutions of the differential equation:  ( 1- x2 ) d2y/dx2 - 2.x.dy/dx + [ n(n+1) - m2/( 1- x2 ) ] y = 0
-Thus, we can seek these solutions by substituting the series  y = a0 + a1.x + a2.x2 + ...... + ak.xk + ...... in this differential equation.
-This approach leads to the recurrence relation:  (k+1)(k+2) ak+2 = [ m2+2k2-n(n+1) ] ak + [ (k-2)(1-k) + n(n+1) ] ak-2        ( a-1= a-2 = 0 )

and  a0 = 2m/pi1/2 cos pi(n+m)/2 . Gam((n+m+1)/2) / Gam((n-m+2)/2)
        a1 = 2m+1/pi1/2 sin pi(n+m)/2 . Gam((n+m+2)/2) / Gam((n-m+1)/2)    for  Pnm(x)

        a0 = -2m-1pi1/2 sin pi(n+m)/2 Gam((n+m+1)/2) / Gam((n-m+2)/2)
        a1 = 2m pi1/2 cos pi(n+m)/2 Gam((n+m+2)/2) / Gam((n-m+1)/2)   for  Qnm(x)
 

Data Registers:      R00 to R08: temp  ( when the program stops, R00 = x & R01 = Pnm(x)  or  Qnm(x)
Flags: F10 & F02    if flag F02 is clear "ALF4" calculates  Pnm(x)
                                 if flag F02 is set    "ALF4" ----------  Qnm(x)
Subroutine:  "GAM"  ( Gamma Function )
 

001  LBL "ALF4"
002  STO 00
003  RDN
004  STO 03
005  X<>Y
006  STO 02
007  +
008  STO 08
009  1
010  +
011  2
012  /
013  STO 04
014  XEQ "GAM"
015  STO 06
016  RCL 03
017  RCL 02
018  -
019  2
020  ST+ Y
021  /
022  STO 05
023  FRC
024  X#0?
025  GTO 02
026  LASTX
027  X>0?
028  GTO 02
029  CLX
030  STO 06
031  GTO 03
032  LBL 02
033  LASTX
034  XEQ "GAM"
035  ST/ 06
036  LBL 03
037  RCL 04
038  .5
039  ST- 05
040  +
041  XEQ "GAM"
042  STO 07
043  RCL 05
044  FRC
045  X#0?
046  GTO 02
047  LASTX
048  X>0?
049  GTO 02
050  CLX
051  STO 07
052  GTO 03
053  LBL 02
054  LASTX
055  XEQ "GAM"
056  ST/ 07
057  LBL 03
058  RCL 08
059  1
060  ASIN
061  *
062  1
063  P-R
064  FS? 02
065  X<>Y
066  FS? 02
067  CHS
068  ST* 06
069  X<>Y
070  ST* 07
071  2
072  RCL 02
073  ST* 02
074  Y^X
075  FC? 02
076  ST+ X
077  PI
078  SQRT
079  FC? 02
080  1/X
081  *
082  ST* 07
083  2
084  /
085  ST* 06
086  RCL 00
087  STO 10
088  RCL 07
089  *
090  RCL 06
091  +
092  STO 01
093  CLX
094  STO 04
095  STO 05
096  STO 08
097  RCL 03
098  1
099  +
100  ST* 03
101  LBL 01
102  3
103  STO 09
104  LBL 04
105  RCL 08
106  2
107  -
108  1
109  RCL 08
110  -
111  *
112  RCL 03
113  +
114  RCL 04
115  *
116  RCL 02
117  RCL 08
118  X^2
119  ST+ X
120  +
121  RCL 03
122  -
123  RCL 06
124  *
125  +
126  RCL 08
127  1
128  ST+ 08
129  +
130  RCL 08
131  LASTX
132  +
133  *
134  /
135  X<> 07
136  X<> 06
137  X<> 05
138  STO 04
139  RCL 00
140  RCL 10
141  *
142  STO 10
143  RCL 07
144  *
145  RCL 01
146  +
147  STO 01
148  LASTX
149  X#Y?
150  GTO 01
151  DSE 09
152  GTO 04
153  END

( 208 bytes / SIZE 010 )
 
 
 
      STACK        INPUTS      OUTPUTS
           Z             m             /
           Y             n             /
           X             x          F(x)

with  F(x) =  Pnm(x)  if flag F02 is clear  and   F(x) =   Qnm(x) if flag F02 is set.
 

Example:    Calculate   P1.30.4(0.7)  &   Q1.30.4(0.7)

 CF 02    0.4   ENTER^
               1.3   ENTER^
               0.7   XEQ "ALF4"  >>>>   P1.30.4(0.7)  =  0.274932821     ( in 1mn56s )

 SF 02    0.4   ENTER^
               1.3   ENTER^
               0.7   XEQ "ALF4"  >>>>   Q1.30.4(0.7)  =  -1.317935686    ( in  1mn43s )

-Lines 23 to 33 ; 36 ; 44 to 54 ; 57 are only useful if    Gam((n-m+2)/2)  or   Gam((n-m+1)/2)   are infinite  ( and the other Gammas are not )
-Thus, we can compute  Q13(0.2) = -1.701034547
-If  (n-m+2)/2  and  (n-m+1)/2  #  0 ; -1 ; -2 ; -3 ; ....  these 24 lines may be deleted.
-An alternative is to use a routine which computes 1/Gam(x) instead of Gam(x): this function is defined for any real x-value!
 

        g) Associated Legendre Functions of the second kind ( fractional order & degree ) with  | x | < 1
 

Formula:      Qnm(x) = 2m pi1/2  (1-x2)-m/2 [ -Gam((1+m+n)/2)/(2.Gam((2-m+n)/2)) . sin pi(m+n)/2 .  F(-n/2-m/2 ; 1/2+n/2-m/2 ; 1/2 ; x2 )
                                                                     + x.Gam((2+n+m)/2) / Gam((1+n-m)/2) . cos pi(m+n)/2 . F((1-m-n)/2 ; (2+n-m)/2 ; 3/2 ; x2 )
 
 

Data Registers:         R00 to R07: temp
Flags: /
Subroutines:  "1/G"  ( 1/Gamma function )
                        "HGF"   ( Hypergeometric functions )
 

01  LBL "ALF5"
02  STO 04
03  RDN
04  STO 05
05  X<>Y
06  STO 06
07  +
08  .5
09  STO 03
10  *
11  CHS
12  STO 01
13  1
14  RCL 05
15  +
16  RCL 06
17  -
18  RCL 03
19  *
20  STO 02
21  RCL 04
22  X^2
23  XEQ "HGF"
24  STO 07
25  RCL 02
26  RCL 03
27  +
28  XEQ "1/G"
29  2
30  /
31  ST* 07
32  RCL 02
33  RCL 06
34  +
35  XEQ "1/G"
36  ST/ 07
37  RCL 03
38  ST+ 01
39  ST+ 02
40  ISG 03
41  RCL 00
42  XEQ "HGF"
43  ST* 04
44  RCL 01
45  RCL 05
46  +
47  XEQ "1/G"
48  ST* 04
49  RCL 02
50  RCL 06
51  +
52  XEQ "1/G"
53  ST/ 04
54  RCL 05
55  RCL 06
56  +
57  1
58  ASIN
59  *
60  SIN
61  ST* 07
62  LASTX
63  COS
64  RCL 04
65  *
66  RCL 07
67  -
68  PI
69  SQRT
70  *
71  4
72  1
73  RCL 00
74  -
75  /
76  RCL 06
77  2
78  /
79  Y^X
80  *
81  END
 

( 125 bytes / SIZE 008 )
 
 
 
      STACK        INPUTS      OUTPUTS
           Z             m             /
           Y             n             /
           X             x         Qnm(x) 

Examples:       0.4  ENTER^   1.3   ENTER^   0.7   XEQ "ALF5"   >>>>     Q1.30.4(0.7) =  -1.317935680   ( in 66 seconds )
                         2.4  ENTER^   0.4   ENTER^   0.1        R/S            >>>>      Q0.42.4(0.1) =   0.102528105   ( in 30 seconds )
 

7°) 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  ""                TEXT0   or another NOP instruction like  LBL 02  STO X   ABS  ...  etc  ...
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

8°) 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  ""                TEXT0   or another NOP instruction like  LBL 02  STO X   ABS  ...  etc  ...
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
 
 

9°) 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 ephemeride 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 )
 

10°) Other Orthogonal Polynomials
 

        a) Ultraspherical Polynomials

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 - ( The modified version above )  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
 

        b) Generalized Laguerre's 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

Note:   Ln(0)(x) = (1/n!).Ln(x)   where   Ln(x) = Laguerre's Polynomial
 

        c) Jacobi's Polynomials
 
 

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
 

11°) Struve Functions
 

Definition:   Hn(x) =  (x/2)n+1 SUM k>=0  (-1)k(x/2)2k/ ( Gam(k+3/2).Gam(k+n+3/2) )     ( here Hn(x) are not  Hermite Polynomials )
 

Data Registers:  /
Flags: /
Subroutine:  "GAM"  ( Gamma function )
 

01  LBL "HNX"
02  2
03  /
04  STO M
05  X<>Y
06  STO N
07  STO O
08  1
09  +
10  Y^X
11  ST+ X
12  PI
13  SQRT
14  /
15  1.5
16  ST+ N
17  ST+ O
18  RCL Y
19  ENTER^
20  LBL 01
21  X<> T
22  RCL M
23  X^2
24  *
25  CHS
26  R^
27  /
28  RCL N
29  /
30  STO T
31  ISG Z
32  ""                TEXT0   or another NOP instruction like  LBL 02  STO X  ...  etc  ...
33  ISG N
34  ""                TEXT0   or another NOP instruction like  LBL 02  STO X  ...  etc  ...
35  X<>Y
36  ST+ Y
37  X#Y?
38  GTO 01
39  X<> O
40  XEQ "GAM"
41  ST/ O
42  RCL O
43  CLA
44  END

( 76 bytes / SIZE 000 )
 
 
      STACK        INPUTS     OUTPUTS
           Y             n    Gam(n+3/2)
           X             x         Hn(x)

Example:    Compute  H1.2 ( 3.4 )

  1.2   ENTER^
  3.4   XEQ "HNX"  yields  1.113372654    ( in 14 seconds )
 

12°) Theta Functions

Formulae:   with  q =  e-pi K'/K  ( see Jacobian elliptic functions for the HP-41 )

   Theta1(x;q) =  2.q1/4 SUMk>=0  (-1)k qk(k+1) sin(2k+1)x
   Theta2(x;q) =  2.q1/4 SUMk>=0    qk(k+1) cos(2k+1)x
   Theta3(x;q) =  1 + 2 SUMk>=1    qk*k cos 2kx                                  ( 0 <= q < 1 )
   Theta4(x;q) =  1 + 2 SUMk>=1  (-1)k qk*k cos 2kx
 

Data Registers:  /
Flags: F01 , F02 , F03 , F04   ( Set  only one of these 4 Flags )
Subroutines: /

-If flag F01 is set , "THETA" calculates Theta1(x;q)
-If flag F02 is set , "THETA" calculates Theta2(x;q)
-If flag F03 is set , "THETA" calculates Theta3(x;q)
-If flag F04 is set , "THETA" calculates Theta4(x;q)

-Synthetic registers M , N , O may be replaced by any unused data registers.
 

01  LBL "THETA"
02  RAD
03  STO M
04  10
05  CHS
06  RCL Z
07  ABS
08  STO N
09  X#0?
10  LOG
11  X#0?
12  /
13  SQRT
14  INT
15  STO O
16  CLX
17  ISG O
18  LBL 01
19  RCL N
20  RCL O
21  FC? 01
22  FS?02
23  ISG X
24  ""                     TEXT0   or another NOP instruction like  LBL 02   STO X   ABS  ...  etc  ...
25  Y^X
26  FC? 01
27  FS? 04
28  CHS
29  RCL O
30  Y^X
31  LASTX
32  ST+ X
33  FC? 01
34  FS? 02
35  ISG X
36  ""                     TEXT0   or another NOP instruction like  LBL 02   STO X   ABS  ...  etc  ...
37  RCL M
38  *
39  FS? 01
40  SIN
41  FC? 01
42  COS
43  *
44  +
45  DSE O
46  GTO 01
47  FC? 01
48  FS? 02
49  RCL M
50  FS? 01
51  SIN
52  FS? 02
53  COS
54  FC? 03
55  FS? 04
56  .5
57  +
58  ST+ X
59  RCL N
60  SQRT
61  SQRT
62  FC? 03
63  FS? 04
64  SIGN
65  *
66  RCL M
67  SIGN
68  X<> N
69  X<>Y
70  DEG
71  CLA
72  END

( 119 bytes / SIZE 000 )
 
 
      STACK        INPUTS     OUTPUTS
           Y             q            q
           X             x    Thetan (x;q)
           L             /            x

Example:   Compute  Theta1(x;q) , Theta2(x;q) , Theta3(x;q) , Theta4(x;q)   for  x = 2 ; q = 0.3

- CF02 CF 03 CF 04

- SF 01                0.3  ENTER^   2   XEQ "THETA"  gives      1.382545289        ( in 12 seconds )
- CF01   SF 02    0.3  ENTER^   2         R/S              -----     -0.488962527         --------------
- CF 02  SF 03    0.3  ENTER^   2         R/S              -----      0.605489938         --------------
- CF 03  SF 04    0.3  ENTER^   2         R/S              -----      1.389795845         --------------
 

13°) Riemann Zeta Function
 

-The Riemann Zeta function is defined by the series:       Zeta(x) = 1 + 1/2x + 1/3x + ........ + 1/nx + ......     ( x>1)
  but the following program uses the formula:                 ( 1 - 1/2x + 1/3x - 1/4x + ...... ) / ( 1 - 21-x )   which converges more rapidly
 

Data Registers:  R00 & R01 ( temp )
Flags: /
Subroutines: /
 

01  LBL "ZETA"
02  CHS
03  STO 01
04  1
05  STO 00
06  ENTER^
07  LBL 01
08  CLX
09  SIGN
10  RCL 00
11  +
12  STO 00
13  RCL 01
14  Y^X
15  -
16  CHS
17  LASTX
18  RND            the accuracy is determined by the display format.
19  X#0?
20  GTO 01
21  SIGN
22  2
23  RCL 01
24  Y^X
25  ST+ X
26  -
27  /
28  ABS
29  END

( 40 bytes / SIZE 002 )
 

Example:    Calculate Apery's number = Zeta(3)

FIX 7     3      XEQ "ZETA"    gives    1.2020569          ( in 3m14s )
FIX 9     3             R/S            -----    1.202056903      ( in 15m12s )

-Execution time and Zeta(x) tend to infinity as x tends to 1.
 

14°) Weierstrass Elliptic Functions  P(x;g2;g3)
 

        a) A Laurent Series ( real arguments )
 

-The Weierstrass Elliptic Function  P(x;g2;g3)  satisfies the differential equation:   (P')2 = 4.P3 - g2.P - g3
-This program calculates  P(x;g2;g3) by a Laurent series:

      P(x;g2;g3) = x-2 + c2.x2 + c3.x4 + ...... + cn.x2n-2 + ....

  where   c2 = g2/20  ;   c3 = g3/28  and   cn =  3 ( c2. cn-2 + c3. cn-3 + ....... + cn-2. c2 ) / (( 2n+1 )( n-3 ))      ( n > 3 )

-The successive  cn  are computed and stored into registers  R05  R06  ......

Data Registers:  R00 to R04 are used for temporary data storage and  R05 = c2 ; R06 = c3 ; R07 = c4 ; ......
Flags: /
Subroutines: /
 

01  LBL "WEF"
02  X^2
03  STO 02
04  CLX
05  20
06  /
07  STO 05
08  X<>Y
09  28
10  /
11  STO 06
12  RCL 02
13  *
14  +
15  STO 01
16  5
17  STO 00
18  LBL 01
19  2
20  STO 04
21  LBL 02
22  RCL 00
23  .1
24  %
25  5
26  +
27  STO 03
28  CLX
29  LBL 03
30  RCL IND Y
31  RCL IND 03
32  *
33  +
34  DSE Y
35  ISG 03
36  GTO 03
37  3
38  *
39  RCL 00
40  ENTER^
41  ST+ Y
42  DSE Y
43  4
44  -
45  *
46  /
47  ISG 03
48  CLX
49  STO IND 03
50  RCL 02
51  RCL 00
52  3
53  -
54  Y^X
55  *
56  RCL 01
57  +
58  STO 01
59  ISG 00
60  CLX
61  LASTX
62  X#Y?
63  GTO 01
64  DSE 04
65  GTO 02
66  RCL 02
67  ST* Y
68  1/X
69  +
70  END

( 95 bytes / SIZE ? )
 
 
 
      STACK        INPUTS     OUTPUTS
           Z             g3            /
           Y             g2            /
           X             x     P(x;g2;g3)

Example:   Calculate   P(x;g2;g3)  for  x = 0.6   g2 = 0.9  g3 = 1.4

  1.4  ENTER^
  0.9  ENTER^
  0.6  XEQ "WEF"  >>>>  2.800500784   ( in 32 seconds )

-The SIZE must be ( at least ) 016 in this example.
-If you get the message "NON EXISTENT" ( line 49 ) , increase the SIZE and R/S
-The Laurent series may diverge if x is too great. In this case, one may use the duplication formula below.
 

        b) Duplication formula ( real arguments )
 

-We have    P(2x) = -2 P(x) + ( 6 P2(x) - g2/2 )2 / ( 4 ( 4 P3(x) - g2 P(x) - g3 ) )
-This formula is used recursively ( n times ) to obtain  P(2nx) if we know  P(x)     ( n = 1 ; 2 ; 3 ; ..... )

Data Registers:   R00 = n  ( n > 0 )   this register is to be initialized before executing "DUPW"
                                R01: temp  R02 = g2  R03 = g3
Flags: /
Subroutine:  /
 

01  LBL "DUPW"
02  RDN
03  STO 02
04  X<>Y
05  STO 03
06  R^
07  LBL 01
08  STO 01
09  X^2
10  6
11  *
12  RCL 02
13  2
14  /
15  -
16  X^2
17  RCL 01
18  ST+ X
19  X^2
20  RCL 02
21  -
22  RCL 01
23  *
24  RCL 03
25  -
26  4
27  *
28  /
29  RCL 01
30  ST+ X
31  -
32  DSE 00
33  GTO 01
34  END

( 47 bytes / SIZE 004 )
 
 
 
      STACK        INPUTS     OUTPUTS
           Z             g3            /
           Y             g2            /
           X      P(x;g2;g3)     P(2nx;g2;g3)

Example:   Calculate   P(x;g2;g3)  for  x = 4.8   g2 = 0.9  g3 = 1.4

-The Laurent series diverges for these arguments but  4.8 = 23* 0.6  and we have already found  that P(0.6;0.9;1.4) = 2.800500784 , thus:

    3   STO 00
  1.4  ENTER^
  0.9  ENTER^
  2.800500784  XEQ "DUPW"  >>>>  P(4.8;0.9;1.4) = 1.954704160  ( in 3seconds )
 

        c) A Laurent Series ( complex arguments )
 

-The same Laurent series is now applied to  z = x + i.y   ( but  g2 and g3 real )

Data Registers:  R00 to R06 are used for temporary data storage and  R07 = c2 ; R08 = c3 ; R09 = c4 ; ......
Flags: /
Subroutines: /
 

001  LBL "WEFZ"
002  R-P
003  X^2
004  STO 03
005  RDN
006  ST+ X
007  STO 04
008  CLX
009  20
010  /
011  STO 07
012  CLX
013  28
014  /
015  STO 08
016  RCL 03
017  *
018  RCL 04
019  X<>Y
020  P-R
021  RCL 07
022  +
023  STO 01
024  X<>Y
025  STO 02
026  7
027  STO 00
028  LBL 01
029  2
030  STO 06
031  LBL 02
032  RCL 00
033  .1
034  %
035  7
036  +
037  STO 05
038  CLX
039  LBL 03
040  RCL IND Y
041  RCL IND 05
042  *
043  +
044  DSE Y
045  ISG 05
046  GTO 03
047  3
048  *
049  RCL 00
050  ENTER^
051  ST+ Y
052  DSE X
053  5
054  ST- Z
055  -
056  *
057  /
058  ISG 05
059  CLX
060  STO IND 05
061  RCL 04
062  RCL 03
063  RCL 00
064  5
065  -
066  ST* Z
067  Y^X
068  RCL IND 05
069  *
070  P-R
071  RCL 01
072  +
073  STO 01
074  LASTX
075  -
076  X<>Y
077  RCL 02
078  +
079  STO 02
080  LASTX
081  -
082  R-P
083  ISG 00
084  CLX
085  X#0?
086  GTO 01
087  RCL 02
088  RCL 01
089  R-P
090  RCL 03
091  *
092  X<>Y
093  RCL 04
094  +
095  X<>Y
096  P-R
097  RCL 04
098  CHS
099  RCL 03
100  1/X
101  P-R
102  RDN
103  ST+ Z
104  X<> T
105  +
106  END

( 135 bytes / SIZE? )
 
 
      STACK        INPUTS     OUTPUTS
           T             g3            /
           Z             g2            /
           Y             y           B
           X             x           A

    with  P(x+iy;g2;g3) = A + i.B

Example:   Calculate   P(z;g2;g3)  for  z = 0.6 + 0.4 i    g2 = 0.9  g3 = 1.4

  1.4  ENTER^
  0.9  ENTER^
  0.4  ENTER^
  0.6  XEQ "WEFZ"  >>>>  0.7390436447    X<>Y    -1.744031391    ( in 48seconds ;  SIZE 018 ( at least ) )

Whence    P(0.6+0.4i ;0.9;1.4) =  0.7390436447 -1.744031391 i
 

        d) Duplication formula ( complex arguments )
 

-The same duplication formula holds for complex arguments.
 

Data Registers:   R00 = n  ( n > 0 )   this register is to be initialized before executing "DUPWZ"
                                R01: temp  R02 = g2  R03 = g3   R04 to R06: temp
Flags: /
Subroutines:  /
 

01  LBL "DUPWZ"
02  R^
03  STO 03
04  R^
05  STO 02
06  R^
07  R^
08  LBL 01
09  STO 01
10  X<>Y
11  STO 04
12  X<>Y
13  R-P
14  3
15  ST* Z
16  Y^X
17  4
18  *
19  P-R
20  RCL 04
21  RCL 02
22  *
23  ST- Z
24  CLX
25  RCL 01
26  LASTX
27  *
28  -
29  RCL 03
30  -
31  R-P
32  STO 05
33  X<>Y
34  STO 06
35  RCL 04
36  RCL 01
37  R-P
38  2
39  ST* Z
40  Y^X
41  6
42  *
43  P-R
44  RCL 02
45  2
46  /
47  -
48  R-P
49  2
50  ST/ Y
51  ST* Z
52  Y^X
53  RCL 06
54  ST- Z
55  CLX
56  RCL 05
57  /
58  P-R
59  RCL 04
60  ST+ X
61  ST- Z
62  CLX
63  RCL 01
64  ST+ X
65  -
66  DSE 00
67  GTO 01
68  END

( 89 bytes / SIZE 007 )
 
 
 
      STACK        INPUTS     OUTPUTS
           T             g3            /
           Z             g2            /
           Y             B           B'
           X             A           A'

    with  P(z;g2;g3) = A + i.B  and   P(2nz;g2;g3) = A' + i.B'
 

Example:     Deduce  P(4.8+3.2 i ;0.9;1.4) from  P(0.6+0.4i ;0.9;1.4) =  0.7390436447 -1.744031391 i
 

     3   STO 00
   1.4  ENTER^
   0.9  ENTER^
   0.4  ENTER^
   0.6  XEQ "DUPWZ"  >>>>  0.152165025   X<>Y   -1.254979593    ( in 22 seconds )

whence   P(4.8+3.2 i ;0.9;1.4) =   0.152165025 - 1.254979593 i
 

        e) Weierstrass & Jacobian Elliptic Functions - Half Periods
 

-The Weierstrass Elliptic Functions may also be computed by mean of the Jacobian Elliptic Functions.
-This program also gives the first derivative and the half-periods of  P(x;g2;g3)
-We have    P(x+2k*Omega+2ik'*Omega') = P(x)   for  all x   ( except the poles )   if  k and k' are integers & Omega and  i.Omega' are the half-periods.

Formulae:      Let  r1 ;  r2 ;  r3   the 3 roots of the polynomial  4x3- g2.x -g3    and   m  the parameter of the complete elliptic integrals of the 1st kind:  K(m)
 

      IF  THE  3 ROOTS  ARE  REAL  (  R1 > R2 > R3 )          IF ONLY ONE ROOT  IS  REAL  (  SAY  R1  )

                              m  =  ( r2 - r3 )/( r1 - r3 )                                                                                           m = 1/2 - 3 r1/(4H)
                            Omega =  K(m)/( r1 - r3 )1/2                                                                                     Omega2 = K(m)/H1/2
                          Omega' =  K(1-m)/( r1 - r3 )1/2                                                                                 Omega'2 = K(1-m)/H1/2

            P(x;g2;g3) = r3 + ( r1 - r3 )/sn2( y | m )                                        P(x;g2;g3) = r1 + H ( 1 + cn(y'|m) ) / ( 1 - cn(y'|m) )
           P'(x;g2;g3) = -2 ( r1 - r3 )3/2 cn (y|m) dn(y|m) / sn3(y|m)                  P'(x;g2;g3) = -4 H3/2 sn(y'|m) dn(y'|m) / ( 1 - cn(y'|m) )2

                           with   y = ( r1 - r3 )1/2.x                                                                           with   y' = 2.x.H1/2    and   H = ( 2.r12 + r2.r3 )1/2
                                                                                                                                                    Omega2 = Omega + i. Omega'
                                                                                                                                                   Omega'2 =  Omega' + i. Omega
 

Data Registers:  R00 to R07 are used for temporary data storage by the Jacobian Elliptic Function and Complete Elliptic Integral programs
    and when the program stops:  R08 = x  ;  R09 = Omega  if  flag F01 is clear  R09 = Omega2  if flag F01 is set      R11: temp
                                                                   R10 = Omega2 -----------------   R10 = Omega'2  ---------------       R12: temp
Flags: F01 - F02 - F05

Subroutines:  "CEI"  ( Complete Ellptic Integrals ) see "Jacobian Elliptic Functions for the HP-41"
                         "JEF" ( Jacobian Elliptic Functions )  ---------------------------------------------
                          "P3"  ( Cubic Equations ) see "Polynomials for the HP-41" or another program which returns the larger real solution of a cubic in X-register.

Note:   If you have used registers R08 & R09 instead of synthetic registers M & N in the "JEF" program,
           replace R08 and R09 by R13 and R14 in the following listing.
 

001  LBL "WEF2"
002  STO 08
003  RDN
004  CHS
005  0
006  X<> Z
007  CHS
008  4
009  RDN
010  XEQ "P3"
011  FS? 01
012  GTO 02
013  RCL Z
014  STO 11
015  ST- Z
016  -
017  STO 12
018  XEQ 01
019  ST/ Y
020  X^2
021  STO 01
022  /
023  *
024  RCL 01
025  1/X
026  GTO 03
027  LBL 01
028  CF 05
029  /
030  STO 10
031  1
032  X<>Y
033  X=Y?
034  SF 05
035  FC? 05
036  XEQ "CEI"
037  X<>Y
038  FS?C 05
039   E90
040  STO 09
041  1
042  STO Y
043  RCL 10
044  -
045  X=Y?
046  SF 05
047  FC? 05
048  XEQ "CEI"
049  X<>Y
050  FS?C 05
051   E90
052  X<> 10
053  RCL 12
054  SQRT
055  ST/ 09
056  ST/ 10
057  RCL 08
058  FS? 01
059  ST+ X
060  *
061  XEQ "JEF"
062  RTN
063  LBL 02
064  STO 11
065  X^2
066  ST+ X
067  RDN
068  R-P
069  X^2
070  R^
071  +
072  SQRT
073  STO 12
074  RCL 11
075  3
076  *
077  X<>Y
078  /
079  2
080  -
081  CHS
082  4
083  XEQ 01
084  SF 01
085  X<>Y
086  STO 01
087  1
088  -
089  STO 02
090  X^2
091  /
092  *
093  ST+ X
094  1
095  RCL 01
096  +
097  RCL 02
098  CHS
099  /
100  LBL 03
101  RCL 12
102  SQRT
103  ST+ X
104  CHS
105  ST* Z
106  X<> L
107  ST* Z
108  *
109  RCL 11
110  +
111  END

( 170 bytes / SIZE 013 )
 
 
      STACK        INPUTS     OUTPUTS
           Z             g3            /
           Y             g2      P'(x;g2;g3)
           X             x      P(x;g2;g3)

Example1:   Calculate   P(x;g2;g3) &  P'(x;g2;g3)  for  x = 2   g2 = 4  g3 = 1

   1   ENTER^
   4   ENTER^
   2  XEQ "WEF2"  >>>>   P(2;4;1) =  4.950267724       X<>Y P'(2;4;1) =  21.55057197    ( in 28 seconds )

-We have   R09 = 1.225694692  &  R10 = 1.496729323   ( Omega & Omega' because F01 is clear )

       Therefore the primitive half-periods are:   1.225694692  &  1.496729323 i
 

Example2:   Calculate   P(x;g2;g3) &  P'(x;g2;g3)  for  x = 1   g2 = 2  g3 = 3

   3   ENTER^
   2   ENTER^
   1  XEQ "WEF2"  >>>>   P(1;2;3) =  1.214433709       X<>Y P'(1;2;3) =  -1.317406193    ( in 27 seconds )

-We have   R09 = 1.197220889  &  R10 = 2.350281226   ( Omega2 & Omega'2 because F01 is set )

       Whence:     Omega =  0.598610445 - 1.175140613 i   &   Omega' =  0.598610445 + 1.175140613 i
 

Notes:

-I've called  Omega' & Omega'2 what Abramowitz and Stegun call    Omega'/i  &  (Omega'2)/i
-This program doesn't work if  g2 = g3 = 0  but in this case,  P(x) = 1/x2
-Lines 28-32 -33-34-35-38-39-42-45-46-47-50-51  are useful only if 2 roots of the polynomial  4x3- g2.x -g3 are equal:
  in this case, one of the half-periods is infinite and a number near E90 is stored in register R09 or R10.
  Without these lines, there would be a DATA ERROR at line 38 of the "CEI" program  ( K(1) = infinite )

  For instance  P(1;48;-64) = 2.181597704 ; P'(1;48;-64) = -0.903006185 ; R09 = 4.082 1089 = infinite = Omega ; R10 = 0.641274915 = Omega'  ( CF01 )
           P(1;12;8) = 2.079381536 ; P'(1;12;8) = 1.735214810 ; R09 = 0.906899682 = Omega ; R10 = 5.7735 1089 = infinite = Omega'  ( CF01 )

-Numerous other relations with Jacobian Elliptic Functions and Theta Functions may be found in the reference below
 

15°) Exponential, Sine, Cosine Integrals

       a) Series Expansions

-These functions are defined by:

          Ei(x) = §-infinityx  et/t dt      ;      Si(x) = §0x sin(t)/t dt      ;      Ci(x) = C + Ln(x) + §0x (cos(t)-1)/t dt

          Shi(x) = §0x sinh(t)/t dt      ;      Chi(x) = C + Ln(x) + §0x (cosh(t)-1)/t dt

-The program computes  Ei(x) , Si(x) , Ci(x) , Shi(x) and Chi(x) with the following formulae:
 

   Ei(x)  = C + ln(x) + Sumn=1,2,.....  xn/(n.n!)
   Si(x)  = Sumn=0,1,2,..... (-1)n x2n+1/((2n+1).(2n+1)!)
  Ci(x)  = C + ln(x) + Sumn=1,2,..... (-1)n x2n/(2n.(2n)!)                 where  C = 0.5772156649 is Euler's constant.
  Shi(x) = Sumn=0,1,2,.....  x2n+1/((2n+1).(2n+1)!)
  Chi(x)= C + ln(x) + Sumn=1,2,.....  x2n/(2n.(2n)!)
 

Data Registers:  R00 = x , R01 to R04: temp
Flags: /
Subroutines: /
 

01  LBL "EI"
02  STO 00
03  STO 04
04  1
05  STO 01
06  STO 03
07  X<>Y
08  GTO 03
09  LBL "SI"
10  STO 00
11  STO 02
12  X^2
13  CHS
14  GTO 01
15  LBL "CI"
16  STO 00
17  X^2
18  CHS
19  GTO 02
20  LBL "SHI"
21  STO 00
22  STO 02
23  X^2
24  LBL 01
25  STO 04
26  1
27  STO 01
28  2
29  STO 03
30  LASTX
31  GTO 04
32  LBL "CHI"
33  STO 00
34  X^2
35  LBL 02
36  STO 04
37  2
38  STO 01
39  STO 03
40  /
41  LBL 03
42  STO 02
43  RCL 01
44  /
45  XEQ 04
46  RCL 00
47  LN
48  .5772156649
49  +
50  +
51  RTN
52  LBL 04
53  RCL 02
54  RCL 04
55  *
56  RCL 01
57  RCL 03
58  ST+ 01
59  SIGN
60  +
61  RCL 01
62  X=Y?
63  SIGN
64  *
65  /
66  STO 02
67  RCL 01
68  /
69  +
70  X#Y?
71  GTO 04
72  END

( 119 bytes / SIZE 005 )
 

Example:    1.4   XEQ "EI"     >>>>  Ei(1.4) =  3.007207463    ( in 8 seconds )
                    1.4   XEQ "SI"     >>>>  Si(1.4) =  1.256226733    ( in 4 seconds )
                    1.4   XEQ "CI"    >>>>  Ci(1.4) =  0.462006585    ( in 5 seconds )
                    1.4   XEQ "SHI"  >>>> Shi(1.4) =  1.561713387   ( in 4 seconds )
                    1.4   XEQ "CHI" >>>> Chi(1.4) =  1.445494076   ( in 5 seconds )

Note:   Logarithmic Integral = Li(x) = Ei(Ln(x))    ( x > 1 )    e.g.   Li(100) = 30.12614160
            ( simply add  LBL "LI"   LN   before the first line if needed )
 

       b) Asymptotic Series
 

-For large arguments, Si(x) and Ci(x) are not calculated with a good accuracy by the above program.
-The following one uses these formulae:

   Si(x) = pi/2 - f(x) cos x - g(x) sin x                  where    f(x) ~  ( 1 - 2!/x2 + 4!/x4 - 6!/x6 + ..... )/x
   Ci(x) = f(x) sin x - g(x) cos x                             and     g(x) ~  ( 1 - 3!/x2 + 5!/x4 - 7!/x6 + ..... )/x2

-Actually, these series are divergent but the error is small if we stop at a proper instant.
 

Data Registers:  R00 = x , R01 to R03: temp
Flags: /
Subroutines: /
 

01  LBL "SCI"
02  STO 00
03  CLX
04  STO 01
05  SIGN
06  STO 02
07  XEQ 01
08  RCL 00
09  /
10  STO 03
11  1
12  STO 01
13  STO 02
14  XEQ 01
15  RCL 00
16  STO Z
17  X^2
18  /
19  RAD
20  P-R
21  RCL 00
22  RCL 03
23  P-R
24  ST+ T
25  X<> Z
26  -
27  1
28  ASIN
29  R^
30  -
31  DEG
32  RTN
33  LBL 01
34  RCL 02
35  RCL 01
36  2
37  +
38  STO 01
39  ENTER^
40  DSE X
41  *
42  *
43  RCL 00
44  X^2
45  CHS
46  /
47  STO 02
48  +
49  X#Y?
50  GTO 01
51  END

( 68 bytes / SIZE 004 )
 
 
      STACK        INPUTS      OUTPUTS
           Y             /           Ci(x)
           X             x           Si(x)

Example:    30   XEQ "SCI"   >>>>   Si(30) =   1.566756540     ( in 14 seconds )
                                  X<>Y               Ci(30) =  -0.033032417

Note:   -This program will not work if x is too small:
           -Divergence already appears too soon for x = 29 but if you add  RND  X<>Y  RND  X<>Y  after line 48
           and run the program in FIX 6 ( for instance ) you will have a good approximation:

-With this modification ( and FIX 6 )    20  XEQ "SCI"  yields   Si(20) = 1.54824168  and  Ci(20) = 0.044419864  ( errors are smaller than 10-7 )
 

       c) Exponential Integrals:  En(x)
 

-We have  En(x) = §1+infinity  t -n e -x.t dt        ( x > 0 ;  n = a positive integer )
-This function is computed by a series expansion:

     En(x) = (-x)n-1 ( -Ln x - gamma + Sumk=1,...,n-1 1/k ) / (n-1)!  - Sumk#n-1 (-x)k / (k-n+1) / k!          where  gamma = Euler's Constant = 0.5772156649...
 

Data Registers:     R00 thru R04: temp
Flags: /
Subroutines: /
 

01  LBL "ENX"
02  STO 00
03  X<>Y
04  STO 01
05  STO 04
06  1
07  STO 02
08  STO 03
09  -
10  X#0?
11  1/X
12  GTO 01
13  LBL 00
14  *
15  +
16  LBL 01
17  RCL 02
18  RCL 00
19  CHS
20  *
21  RCL 03
22  /
23  STO 02
24  ISG 03
25  CLX
26  RCL 03
27  RCL 01
28  -
29  X=0?
30  GTO 00
31  /
32  -
33  X#Y?
34  GTO 01
35  0
36  LBL 02
37  DSE 04
38  X<0?
39  GTO 03
40  RCL 04
41  1/X
42  +
43  GTO 02
44  LBL 03
45  RCL 00
46  LN
47  .5772156649
48  +
49  -
50  RCL 00
51  CHS
52  RCL 01
53  DSE X
54  ""                         ( TEXT 0 or another NOP instruction like STO X ... )
55  Y^X
56  LASTX
57  FACT
58  /
59  *
60  +
61  END

( 87 bytes / SIZE 005 )
 
 
      STACK        INPUTS      OUTPUTS
           Y             n             /
           X             x          En(x)

   where  x > 0  &  n is a positive integer

Example:        2    ENTER^
                      1.4  XEQ "ENX"  >>>>  E2(1.4) = 0.083889926      ( in 10 seconds )
 
 

16°) Fresnel Integrals
 

       a) Ascending Series
 

-Fresnel integrals are defined by:    S(x) = §0x  sin(pi.t2/2).dt      and      C(x) = §0x  cos(pi.t2/2).dt

Formulae:       S(x) = SUMn=0,1,2,.....  (-1)n (pi/2)2n+1 x4n+3 / ((2n+1)!(4n+3))

                        C(x) = SUMn=0,1,2,.....  (-1)n (pi/2)2n x4n+1 / ((2n)!(4n+1))
 

Data Registers:  R00 = x ; R01 to R03: temp
Flags: /
Subroutines: /
 

01  LBL "SX"
02  STO 00
03  3
04  Y^X
05  PI
06  6
07  /
08  *
09  1
10  STO 02
11  GTO 00
12  LBL "CX"
13  STO 00
14  0
15  STO 02
16  LBL 00
17  RCL 00
18  X^2
19  PI
20  *
21  2
22  /
23  X^2
24  CHS
25  STO 03
26  R^
27  STO 01
28  LBL 01
29  RCL 01
30  RCL 02
31  2
32  +
33  STO 02
34  ST+ X
35  3
36  -
37  ST* Y
38  4
39  +
40  RCL 02
41  ST* Y
42  DSE X
43  *
44  /
45  RCL 03
46  *
47  STO 01
48  +
49  X#Y?
50  GTO 01
51  END

( 69 bytes / SIZE 004 )
 
 
 
      STACK        INPUTS      OUTPUTS
           X             x      S(x) or C(x)

Example:      2  XEQ "SX"  >>>>  S(2) = 0.343415677    ( in 14 seconds )
                      2  XEQ "CX" >>>>  C(2) = 0.488253412      --------------
 

Note:   This method doesn't give accurate results for large arguments:

              S(3) and C(3) are correct to 5 decimals
           S(3.5) and C(3.5) -----------  4 --------
              S(4) and C(4) are not exact even to 1D!
 

       b) Asymptotic Series
 

Formulae:       S(x) = 1/2 - f(x).cos(pi.x2/2) - g(x).sin(pi.x2/2)     &     C(x) = 1/2 + f(x).sin(pi.x2/2) - g(x).cos(pi.x2/2)

       where               f(x) ~  ( 1 - 1*3/(pi.x2)2 + 1*3*5*7/(pi.x2)4 - ....... ) / (pi.x)
                               g(x) ~  ( 1 - 1*3*5/(pi.x2)2 + 1*3*5*7*9/(pi.x2)4 - ....... ) / (pi2.x3)
 

Data Registers: R00 = x ; R01 to R04: temp
Flags: /
Subroutines: /
 

01  LBL "FSC"
02  STO 00
03  X^2
04  PI
05  *
06  STO 03
07  SIGN
08  STO 01
09  LASTX
10  1/X
11  STO 02
12  XEQ 01
13  STO 04
14  1
15  CHS
16  STO 01
17  CHS
18  STO 02
19  XEQ 01
20  PI
21  RCL 00
22  *
23  ST/ 04
24  /
25  RCL 03
26  2
27  /
28  STO 03
29  X<>Y
30  RAD
31  P-R
32  RCL 03
33  RCL 04
34  P-R
35  ST- T
36  RDN
37  +
38  .5
39  ST+ Z
40  X<>Y
41  -
42  DEG
43  RTN
44  LBL 01
45  RCL 02
46  RCL 01
47  2
48  +
49  ST* Y
50  LASTX
51  +
52  STO 01
53  *
54  RCL 03
55  X^2
56  /
57  CHS
58  STO 02
59  +
60  X#Y?
61  GTO 01
62  END

( 80 bytes / SIZE 005 )
 
 
 
      STACK        INPUTS      OUTPUTS
           Y             /           C(x)
           X             x           S(x)

Examples:        10   XEQ FSC  >>>>   S(10) = 0.468169979    X<>Y    C(10) =  0.499898695    ( in 5 seconds )
                          4.1  XEQ FSC  >>>>   S(4.1) = 0.475798258    X<>Y    C(4.1) = 0.573695632    ( in 14 seconds )
 

-The series are already divergent too soon for  x = 4   but if you add     RND  X<>Y   RND   X<>Y   after line59, you'll get:

         FIX 9     4   XEQ FSC   >>>>   S(4) = 0.420515754   X<>Y   C(4) = 0.498426033   ( in 11 seconds )
 and  FIX 5     3   XEQ FSC   >>>>   S(3) = 0.4963140       X<>Y   S(3) = 0.6057203   ( exact values are  0.4963130  and  0.6057208 to 7D )
 
 

17°) Error Function
 

-This function is defined by:        erf(x) =  (2/pi1/2)  §0x  e-t^2 dt

-To achieve the best accuracy, the following program calculates erf(x) by a series expansion for x < 2 and 1 - erf(x) by continued fractions otherwise  ( x >0 )

Formulae:      erf(x) = (2/pi1/2)  SUMn=0,1,2,.....  (-1)n x2n+1 / (n! (2n+1))

                      2.ex^2  §xinfinity e-t^2 dt = 1/(x+0.5/(x+1/(x+1.5/(x+2/(x+....)))))
 
 

Data Registers: /
Flags: /
Subroutines: /

-Synthetic registers M , N , O  may of course be replaced by any data registers.
 

01  LBL "ERF"
02  ABS
03  STO M
04  STO O
05  2
06  X<=Y?
07  GTO 02
08  SIGN
09  X<>Y
10  ENTER^
11  X^2
12  STO N
13  LBL 01
14  CLX
15  RCL N
16  RCL O
17  *
18  CHS
19  R^
20  ISG T
21  CLX
22  /
23  STO O
24  R^
25  ST+ X
26  DSE X
27  /
28  X<>Y
29  ST+ Y
30  X#Y?
31  GTO 01
32  ST+ X
33  PI
34  SQRT
35  /
36  ENTER^
37  CHS
38  1
39  +
40  X<>Y
41  GTO 04
42  LBL 02
43  X<>Y
44  STO Y
45  24
46  STO N
47  1
48  +
49  X<>Y
50  ST+ X
51  /
52  LBL 03
53  +
54  RCL N
55  X<>Y
56  ST+ X
57  /
58  DSE N
59  GTO 03
60  +
61  X<>Y
62  X^2
63  E^X
64  *
65  1/X
66  PI
67  SQRT
68  /
69  ENTER^
70  CHS
71  1
72  +
73  LBL 04
74  RCL M
75  SIGN
76  RDN
77  CLA
78  END

( 109 bytes / SIZE 000 )
 
 
 
      STACK        INPUTS      OUTPUTS
           Y             /        1-erf(x)
           X          x > 0          erf(x)
           L             /             x

 

Example:         0.9    XEQ ERF  >>>>   erf(0.9) = 0.7969082124      X<>Y     1-erf(0.9) = 0.2030917876            ( in 9 seconds )
                         2.7         R/S       >>>>   erf(2.7) = 0.9998656673      X<>Y     1-erf(2.7) = 0.0001343327399      ( in 7 seconds )
 

18°) Coulomb Wave Functions
 

-The regular Coulomb wave function  FL(n;r)  and the irregular Coulomb wave function GL(n;r) are 2 independant solutions
  of the differential equation:  d2y/dr2 + (1-2n/r-L(L+1)/r2).y = 0
 

       a) Regular Coulomb Wave Function
 

Formulae:         FL(n;r) = CL(n) r L+1 Sum k>L AkL (n) r k-L-1

    with   CL(n) = (1/Gam(2L+2)) 2L e -pi.n/2 | Gam(L+1+i.n) |
    and    AL+1L = 1 ;  AL+2L = n/(L+1)  ;  (k+L)(k-L-1) AkL = 2n Ak-1L - Ak-2L   ( k > L+2 )
 

Data Registers:   R00 thru R05: temp
                              R06 = r  ;  R07 = n  ;  R08 = L  ;  R09 = CL(n)
Flags: /
Subroutine:   "GAMZ"  ( Gamma Function in the complex plane )
 

01  LBL "RCWF"
02  STO 06
03  RDN
04  STO 07
05  X<>Y
06  STO 08
07  1
08  +
09  XEQ "GAMZ"
10  LASTX
11  RCL 07
12  PI
13  *
14  2
15  /
16  E^X
17  /
18  2
19  RCL 08
20  Y^X
21  *
22  RCL 08
23  ST+ X
24  1
25  +
26  FACT
27  /
28  STO 09
29  CLX
30  STO 01
31  SIGN
32  STO 02
33  STO 03
34  LBL 01
35  2
36  STO 04
37  X<>Y
38  LBL 02
39  RCL 07
40  ST+ X
41  RCL 02
42  ST* Y
43  X<> 01
44  RCL 06
45  *
46  -
47  RCL 06
48  *
49  RCL 03
50  ST/ Y
51  RCL 08
52  ST+ X
53  +
54  1
55  ST+ 03
56  +
57  /
58  STO 02
59  +
60  X#Y?
61  GTO 01
62  DSE 04
63  GTO 02
64  RCL 09
65  *
66  RCL 06
67  ST* Y
68  RCL 08
69  Y^X
70  *
71  END

( 96 bytes / SIZE 010 )
 
      STACK        INPUTS      OUTPUTS
           Z             L              /
           Y             n              /
           X             r           FL(n;r)

  where  L is a non-negative integer,
             n  is real,
   and    r  is positive.

Example:

     2   ENTER^
   0.7  ENTER^
   1.8  XEQ "RCWF"  >>>>   F2( 0.7 ; 1.8 ) = 0.141767744   ( in 34 seconds )
 
 

       b) Irregular ( Logarithmic ) Coulomb Wave Function
 
 

Formulae:         GL(n;r) = (1/Pi).(exp(2.Pi.n) - 1) FL(n;r) [ Ln(2r) + qL(n)/pL(n) ] + (1/CL(n)/rL/(2L+1))  Sum k>-L-1  akL (n) r k+L

    with      a -L-1L = 0 ;  a -LL = 1  ;  (k+L)(k-L-1) akL = 2n ak-1L - ak-2L - (2k-1).pL(n).AkL    and   a L+1L = 0

                pL(n) = 2n ( 1+n2 )( 4+n2 ) ........... ( L2+n2 ) 22L / ( 2L+1) / [(2L)!]2

       qL(n)/pL(n) = Sum s=1,...,L s/(s2+n2)  - 1/1 - 1/2 - 1/3 - ...... - 1/(2L+1) + Ré ( Psi(1+i.n) ) + 2.gamma + rL(n)/pL(n)   ( gamma = Euler's Constant )

  where    rL(n) = ( (-1)L+1 / (2L)! )  Im [ 1/(2L+1) + 2(i.n-L)/1!/(2L) + 22(i.n-L)(i.n-L+1)/2!/(2L-1)
                                                               + ...................................... + 22L(i.n-L)(i.n-L+1) ....... (i.n+L-1)/(2L)! ]
 

Data Registers:   R00 thru R05 & R11 thru R15: temp
                              R06 = r  ;  R07 = n  ;  R08 = L  ;  R09 = CL(n)  ;  R10 = FL(n;r)
Flags: /
Subroutines:  "RCWF"  ( regular Coulomb wave function )
                         "PSIZ"    ( Digamma function in the complex plane )
 

  01  LBL "ICWF"
  02  XEQ "RCWF"
  03  STO 10
  04  RCL 07
  05  ST+ X
  06  STO 15
  07  PI
  08  *
  09  E^X-1
  10  *
  11  PI
  12  /
  13  STO 11
  14  RCL 07
  15  1
  16  XEQ "PSIZ"
  17  STO 12
  18  RCL 08
  19  ST+ X
  20  STO 03
  21  1
  22  +
  23  LBL 01
  24  1/X
  25  ST- 12
  26  LASTX
  27  DSE X
  28  GTO 01
  29  STO 02
  30  STO 04
  31  RCL 15
  32  4
  33  RCL 08
  34  Y^X
  35  *
  36  1
  37  STO 13
  38  RCL 03
  39  ST+ Y
  40  FACT
  41  STO 05
  42  X^2
  43  *
  44  /
  45  STO 00
  46  RCL 08
  47  STO 01
  48  X=0?
  49  GTO 00
  50  LBL 02
  51  RCL 01
  52  ENTER^
  53  X^2
  54  RCL 07
  55  X^2
  56  +
  57  ST* 00
  58  /
  59  ST+ 12
  60  DSE 01
  61  GTO 02
  62  LBL 03
  63  RCL 07
  64  RCL 08
  65  RCL 03
  66  -
  67  R-P
  68  X<>Y
  69  RCL 02
  70  +
  71  STO 02
  72  X<>Y
  73  RCL 13
  74  *
  75  ST+ X
  76  RCL 08
  77  ST+ X
  78  RCL 03
  79  -
  80  1
  81  +
  82  /
  83  STO 13
  84  P-R
  85  CLX
  86  RCL 03
  87  /
  88  ST+ 04
  89  DSE 03
  90  GTO 03
  91  LBL 00
  92  RCL 12
  93  1
  94  CHS
  95  RCL 08
  96  Y^X
  97  RCL 04
  98  *
  99  RCL 05
100  /
101  RCL 00
102  X=0?
103  SIGN
104  /
105  -
106  1.15443133       2*Euler's Constant
107  +
108  RCL 06
109  ST+ X
110  LN
111  +
112  ST* 11
113  CLX
114  STO 02
115  STO 03
116  SIGN
117  STO 04
118  STO 05
119  STO 13
120  LBL 04
121  2
122  STO 14
123  LBL 05
124  RCL 15
125  RCL 04
126  ST* Y
127  X<> 03
128  RCL 06
129  *
130  -
131  RCL 06
132  *
133  RCL 15
134  RCL 02
135  ST* Y
136  X<> 01
137  RCL 06
138  *
139  -
140  RCL 06
141  *
142  RCL 05
143  ST/ Y
144  RCL 08
145  ST+ X
146  -
147  1
148  -
149  STO 12
150  X#0?
151  GTO 00
152  RDN
153  CLX
154  RCL 06
155  RCL 05
156  Y^X
157  1
158  LBL 00
159  /
160  STO 02
161  RCL 00
162  *
163  RCL 05
164  RCL 08
165  -
166  ST+ X
167  1
168  -
169  *
170  -
171  RCL 05
172  /
173  RCL 12
174  X#0?
175  /
176  STO 04
177  ISG 05
178  ""                 ( TEXT 0  or another  NOP  instruction like STO X )
179  RCL 13
180  +
181  STO 13
182  LASTX
183  X#Y?
184  GTO 04
185  DSE 14
186  GTO 05
187  RCL 06
188  RCL 08
189  Y^X
190  RCL 09
191  *
192  RCL08
193  ST+ X
194  1
195  +
196  *
197  /
198  RCL 11
199  +
200  END

( 259 bytes / SIZE 016 )
 
 
      STACK        INPUTS      OUTPUTS
           Z             L              /
           Y             n              /
           X             r           GL(n;r)

  where  L is a non-negative integer,
             n  is real,
   and    r  is positive.         FL(n;r) is stored in R10

Example:

            3   ENTER^
         -0.4  ENTER^
          1.2  XEQ "ICWF"  >>>>  G3( -0.4 ; 1.2 ) = 6.563265438        ( in 102 seconds )
                               and we have  F3( -0.4 ; 1.2 ) = 0.029547268          in register R10.

-Unfortunately, several significant digits are sometimes cancelled when the last operation ( line 199 ) is performed.
 ( Check the value in LASTX or in R11 )
 

19°)  Debye Functions
 

-The following program computes   db(x;n) = §x+infinity  tn/(et-1).dt      where n is a positive integer  and  x > 0

Formula:   db(x;n) = Sum k>0  e-k.x [ xn/k + n.xn-1/k2 + ..... + n!/kn+1 ]
 

Data Registers: /
Flags: /
Subroutines: /
 

01  LBL "DEBYE"
02  CLA
03  STO M
04  X<>Y
05  STO N
06  CLST
07  LBL 01
08  R^
09  1
10  +
11  RCL M
12  RCL N
13  STO P           ( synthetic )
14  Y^X
15  RCL Y
16  /
17  ENTER^
18  LBL 02
19  RCL P
20  *
21  R^
22  /
23  RCL M
24  /
25  ST+ Y
26  DSE P
27  GTO 02
28  X<> T
29  RCL M
30  *
31  E^X
32  /
33  RCL O
34  +
35  STO O
36  LASTX
37  X#Y?
38  GTO 01
39  RCL M
40  SIGN
41  X<> N
42  X<>Y
43  CLA
44  END

( 72 bytes / SIZE 000 )
 
 
 
      STACK        INPUTS      OUTPUTS
           Y             n             n
           X             x         db(x,n)
           L             /             x

  n = a positive integer ; x > 0

Example:

    3   ENTER^
  0.7  XEQ "DEBYE"  >>>>  db( 0.7 ; 3 ) = 6.406833597   ( 55 seconds )

Note:    We also have  db(0;n) = §0+infinity  tn/(et-1).dt  =  n!  Zeta(n+1)   where  "Zeta" is the Riemann Zeta Function.
 

20°)  Dawson's Integrals
 

-This program computes   F(x) = e -x^2 §0x  et^2 dt     by a series expansion.
 

Data Registers: /
Flags: /
Subroutines: /
 

01  LBL "DAW"
02  STO M
03  X^2
04  STO N
05  1
06  LASTX
07  ENTER^
08  LBL 01
09  CLX
10  RCL M
11  RCL N
12  *
13  R^
14  ISG T
15  CLX
16  /
17  STO M
18  R^
19  ST+ X
20  DSE X
21  /
22  X<>Y
23  ST+ Y
24  X#Y?
25  GTO 01
26  RCL N
27  E^X
28  /
29  CLA
30  END

( 49 bytes / SIZE 000 )
 
 
 
      STACK        INPUTS      OUTPUTS
           X             x          F(x)
           L             /           e x^2

Examples:

     1.94  XEQ "DAW"  >>>>  F( 1.94 ) = 0.3140571659     ( 11 seconds )
      10         R/S            >>>>    F(10)    = 0.05025384716   ( 85 seconds )
 
 

Reference:   Abramowitz and Stegun , "Handbook of Mathematical Functions" -  Dover Publications -  ISBN  0-486-61272-4
 
 
 

Go back to the HP-41 software library
Go back to the general software library
Go back to the main exhibit hall