The Museum of HP Calculators
The Museum of HP Calculators
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