Special Functions for the HP48S/SX/G/GX

Overview

1°) Gamma Function
2°) Digamma Function
3°) Polygamma Function
4°) Zeta Function
5°) Kummer's Function
6°) Gauss Hypergeometric Functions
7°) Regularized Hypergeometric Functions
8°) Hypergeometric Function U
9°) Generalized Hypergeometric Functions
10°) Appell Hypergeometric Functions

a) Appell Function F₁
b) Appell Function F₂
c) Appell Function F₃
d) Appell Function F₄

11°) Kampé de Fériet Functions
12°) Lauricella Functions

a) Lauricella Function F_A
b) Lauricella Function F_B
c) Lauricella Function F_C
d) Lauricella Function F_D

13°) Srivastava Functions
14°) Associated Legendre Functions
15°) Bessel Functions

a) Bessel & Modified Bessel functions of the 1st & 2nd kind
b) Spherical Bessel Functions
c) Kelvin Functions

16°) Orthogonal Polynomials & Functions

a) Laguerre Polynomials
b) Generalized Laguerre Polynomials & Functions
c) Legendre Polynomials
d) Jacobi Polynomials & Functions
e) Hermite Polynomials & Functions
f) Chebyshev Polynomials & Functions
g) Ultraspherical Polynomials & Functions

17°) Elliptic Functions

a) Weierstrass Elliptic Functions
b) Jacobian Elliptic Functions
c) Carlson Elliptic Integrals
d) Surface Area of an Ellipsoid

18°) Miscellaneous

a) Airy
b) Anger & Weber
c) Dawson Integrals
d) Exponential Integral En
e) Exponential Integral Ei
f) Sine & Cosine Integrals
g) Hyperbolic Sine & Cosine Integrals
h) Catalan Functions
i) Fresnel Integrals
j) Debye Function
k) Parabolic Cylinder Functions
l) Error Function
m) Fibonacci Functions
n) Struve Functions
o) Lommel's Functions
p) PolyLogarithm
q) Incomplete Beta
r) Incomplete Gamma
s) Lerch Transcendent Functions
t) Coulomb Wave Functions
u) Whittaker's Functions
v) Toronto Functions

19°) SpheroidalWave Functions
20°) Ellipsoidal Wave Functions
21°) Euler Gamma Constant
22°) Asymptotic Expansions

a) Airy Functions
b) Coulomb Wave
c) Parabolic Cylinder
d) Exponential Integrals En
e) Struve Function Hn
f) Bessel Functions J & Y & K

23°) Fractional-Integro Differentiation

a) Exponential
b) Logarithm
c) FClog
d) Cosine
e) Sine
f) Hyperbolic Cosine
g) Hyperbolic Sine
h) Airy Functions
i) Error Function
j) Hermite Function
k) Kummer's Function
l) Generalized Laguerre's Function
m) Bessel Functions
n) Fresnel Integrals
o) Exponential Integrals & Related Functions

24°) Non-Linear Equations

a) Secant Method
b) Newton's Method

25°) Continued Fractions

-Most of these programs work with complex arguments

1°) Gamma Function

"GAM" employs Gam(x) = [ x^(x-1/2) ] sqrt(2.Pi) exp [ -x + ( 1/12 )/( x + ( 1/30 )/( x + ( 53/210)/( x + (195/371)/( x + ... )))) ]

and the reccrence relation: Gam(x+n) = (x+n-1) (x+n-2) ...... (x+1) x Gam(x)

STACK	INPUT	OUTPUT
Level 1	x	Gamma(x)

Example:

3.14 GAM -> 2.28448063382

Note:

-This program also works with complex numbers:

(1,2) GAM -> (.151904002671,1.98048801618E-2)

-So Gam ( 1+2.i ) = 0.151904002671 + 0.0198048801618 i

2°) Digamma Function

PSI uses the asymptotic expansion: Psi(x) = ln x - 1/(2x) -1/(12x²) + 1/(120x⁴) - 1/(252x⁶) + 1/(240x⁸) + .......

together with the property: Psi(x+1) = Psi(x) + 1/x

STACK	INPUT	OUTPUT
Level 1	x	Psi(x)

Examples:

3.14 PSI -> .97661709147
(1,2) PSI -> (.71459151538,1.32080728265)

3°) Polygamma Functions

PSIN 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⁽ⁿ⁾(x+1) = Psi⁽ⁿ⁾(x) + (-1)ⁿn! x^-n-1 is used.

STACK	INPUTS	OUTPUTS
Level 2	n	/
Level 1	x	Psi⁽ⁿ⁾(x)

Examples:

1 ENTER
3.14 PSIN -> .374464556478 ( Trigamma )

3 ENTER
(1,2) PSIN -> ( -0.18478233339 , 0.175169315188 ) ( Pentagamma )

4°) Zeta Function

ZETA employs Borwein algorithm ( cf reference [4] ) if Re(z) > 0.5
and the formula: Zeta(z) = Zeta(1-z) Pi^z-1/2 Gamma((1-z)/2) / Gamma(z/2)

STACK	INPUT	OUTPUT
Level 1	x	Zeta(x)

Example:

-7.49 ZETA -> Zeta(-7.49) = 3.31204016904 E-3
(1,2) ZETA -> Zeta( 1+2i ) = 0.598165569763 - 0.351854745218 i

5°) Kummer's Functions

-Kummer's function M(a,b,x) is defined by

M(a;b;x) = 1 + (a)₁/(b)₁. x¹/1! + ............. + (a)_n/(b)_n . xⁿ/n! + .......... where (a)_n = a(a+1)(a+2) ...... (a+n-1)

STACK	INPUTS	OUTPUTS
Level 3	a	/
Level 2	b	/
Level 1	x	M(a,b,x)

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

2
3
PI CHS KUM yields 0.166374561777

KUM also works with complex numbers, for example:

KUM ( 2 , 3 ; 1+2i ) = 0.22245799696 + 1.80488317429 i

Notes:

a)   2x (Pi)^-1/2 M(1/2;3/2;-x²) = erf(x) = error function
b)        (x/2)ⁿ M(n+1/2;2n+1;2x) = Gamma(1+n) e^x I_n(x)        where     I_n = a modified Bessel function
c)       (x^a/a)   M(a;a+1;-x) = incgam(a;x) = §₀^x e^-t t^a-1 dt                ( incgam = incomplete gamma function )

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

-Kummer functions are in fact a special case of the hypergeometric functions

6°) Gauss Hypergeometric Functions

Formulas: F(a,b;c;x) = 1 + (a)₁(b)₁/(c)₁. x¹/1! + ............. + (a)_n(b)_n/(c)_n . xⁿ/n! + .......... where (a)_n = a(a+1)(a+2) ...... (a+n-1) and | x | < 1

-We also have the folowing properties:

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 x < 0 F(a,b,c,x) = ( 1 - x )^-a F(a,c-b,c,-x/(1-x))

If x > 1 F(a,b,c,x) = [ Gam(c) Gam(b-a) / Gam(b) / Gam(c-a) ] (-x)^-a F(a,1+a-c,1+a-b,1/x)
+ [ Gam(c) Gam(a-b) / Gam(a) / Gam(c-b) ] (-x)^-b F(b,1+b-c,1+b-a,1/x) which may be a complex number

STACK	INPUTS	OUTPUTS
Level 4	a	/
Level 3	b	/
Level 2	c	/
Level 1	x	F(a,b;c;x)

Examples: Calculate F(0.3 , -0.7 ; 0.4 ; 0.49)

0.3    ENTER
-0.7   ENTER
   0.4   ENTER
0.49   HGF -> F(0.3 , -0.7 ; 0.4 ; 0.49) = 0.720164355556

-Likewise, we find:

F(1;2;7;1) = 1.5
F(2;3;3;0.7) = 11.1111111111

F(1.2,2.3;3.7;-3) = 0.309850660994
F(1.2,2.3;3.7;4) = -0.804024119422 - 0.105379848361 i

Note:

HGF also works with complex arguments

7°) Regularized Hypergeometric Functions

Formula:

F^~( a , b; c; x ) = SUM_{k=0,1,2,.....} [ (a)_k(b)_k] / [ Gam(k+c) ] . x^k/k!

STACK	INPUTS	OUTPUTS
Level 4	a	/
Level 3	b	/
Level 2	c	/
Level 1	x	F^~( a , b; c; x )

Examples: Calculate F^~(0.3 , -0.7 ; 0.4 ; 0.49)

0.3    ENTER
-0.7   ENTER
   0.4   ENTER
0.49   RHGF -> F^~(0.3 , -0.7 ; 0.4 ; 0.49) = 0.324667518882 it's simply   F(0.3 , -0.7 ; 0.4 ; 0.49) / Gamma ( 0.4 )

-Likewise, F^~(0.3 , -0.7 ; -2 ; 0.49) = -0.0168034752793

8°) Hypergeometric Function U

U(a,b;x) = [ Gam(b-1) / Gam(a) ] x^1-b ₁F₁(a-b+1;2-b;x) + [ Gam(1-b) / Gam(a-b+1) ] ₁F₁(a;b;x) if b is not an integer

U(a;n;x) = [(-1)ⁿ/Gam(a-n+1)] [ Ln(x) / (n-1) ! ₁F₁(a;n;x) + Sum_k=0,1,..... { (a)_k ( Psi(a+k) - Psi(k+1) - Psi(k+n) ) x^k / (k+n-1)! / k! } - Sum_k=1,...,n-1 { (k-1)! x^-k / (1-a)_k /(n-k-1)! } ] if b = n = positive integer

U(a;1-n;x) = xⁿ U(a+n;1+n;x)

U(a,a+n;x) = x^-a SUM_{k=0,.....,n-1} Comb(n-1;n-k-1) (a)_k x^-k if n is a positive integer

STACK	INPUTS	OUTPUTS
Level 3	a	/
Level 2	b	/
Level 1	x	/

Examples:

U(1.2;2.3;3.4) = 0.23681859297

U(1.2;3;1/Pi) = 13.303479155

U(1.2;-3;1/Pi) = 0.163974005823

U(2;7;3.14) = 1.16985623465

9°) Generalized Hypergeometric Functions

-The definition of F(a,b,c;x) may be generalized as follows: given 2 non-negative integers m & p ( at least one of them positive )

_mF_p( a₁,a₂,....,a_m; b₁,b₂,....,b_p; x ) = SUM_{k=0,1,2,.....} [(a₁)_k(a₂)_k.....(a_m)_k] / [(b₁)_k(b₂)_k.....(b_p)_k] . x^k/k!

where (a_i)_k = a_i(a_i+1)(a_i+2) ...... (a_i+k-1) & (a_i)₀ = 1 , likewise for (b_j)_k

-The regularized generalized hypergeometric function F tilde is defined by

_mF^~_p( a₁,a₂,....,a_m; b₁,b₂,....,b_p; x ) = SUM_{k=0,1,2,.....} [ (a₁)_k(a₂)_k.....(a_m)_k] / [Gam(k+b₁) Gam(k+b₂).....Gam(k+b_p)] . x^k/k!

-If no b_j is a negative integer, we simply divide F by the product Gam(b₁) Gam(b₂).....Gam(b_p)
-Otherwise, the calculation is more complicated.

GHGF calculates F or F tilde according to the sign of level2-input.

-All the coefficients must be placed in the stack

STACK	INPUTS	OUTPUTS
...........	a₁	/
...........	..........	/
...........	a_m	/
............	b₁	/
............	.........	/
Level 4	b_p	/
Level 3	m	/
Level 2	+/-p	/
Level 1	x	/

Example1: ₃F₄( 1 , 4 , 7; 2 , 3 , 6 , 5; Pi ) = ?

1 ENTER
4 ENTER
7 ENTER
2 ENTER
3 ENTER
6 ENTER
5 ENTER
3 ENTER
4 ENTER
Pi GHGF -> ₃F₄( 1 , 4 , 7; 2 , 3 , 6 , 5; Pi ) = 1.63101964329

Example2: ₃F^~₄( 1 , 4 , 7; 2 , 3 , 6 , 5; Pi ) = ?

1 ENTER
4 ENTER
7 ENTER
2 ENTER
3 ENTER
6 ENTER
5 ENTER
3 ENTER
-4 ENTER
Pi GHGF -> ₃F^~₄( 1 , 4 , 7; 2 , 3 , 6 , 5; Pi ) = 0.000283163132516

Notes:

-If m = p = 0 GHGF returns exp(x)

-The program doesn't check if the series are convergent or not.
-Even when they are convergent, execution time may be very long !

10°) Appell Functions

a) Appell Function F₁

-The Appell hypergeometric function F₁ is defined by:

F₁( a , b , b' , c ; x , y ) = Sum_{m,n=0 to
infinity} [ (a)_m+n (b)_m (b')_n / (c)_m+n ] ( x^m / m! ) ( yⁿ /n! )

where (a)_k = a(a+1)(a+2).....(a+k-1) if k > 0 and (a)₀ = 1 ( Pochhammer symbol )

STACK	INPUTS	OUTPUTS
Level 6	a	/
Level 5	b	/
Level 4	b'	/
Level 3	c	/
Level 2	x	/
Level 1	y	F₁(a,b,b',c;x,y)

Example:

2 ENTER
3    ENTER
4    ENTER
7 ENTER
   0.1 ENTER
0.2 APL1 -> F₁( 2 , 3 , 4 , 7 , 0.1 , 0.2 ) = 1.40945520741

Note:

-The series is convergent if Max { | x | , | y | } < 1

b) Appell Function F₂

-The Appell hypergeometric function F₂ is defined by:

F₂( a , b , b' , c , c' ; x , y ) = Sum_{m,n=0 to infinity} [ (a)_m+n (b)_m (b')_n / (c)_m (c')_n ] ( x^m / m! ) ( yⁿ /n! )

STACK	INPUTS	OUTPUTS
Level 7	a	/
Level 6	b	/
Level 5	b'	/
Level 4	c	/
Level 3	c'	/
Level 2	x	/
Level 1	y	F₂(a,b,b',c,c';x,y)

Example:

2 ENTER
3    ENTER
4    ENTER
6 ENTER
7 ENTER
   0.1 ENTER
0.2 APL2 -> F₂( 2 , 3 , 4 , 6 , 7 , 0.1 , 0.2 ) = 1.44155353037

Note:

-The series is convergent if | x | + | y | < 1

c) Appell Function F₃

-The Appell hypergeometric function F₃ is defined by:

F₃( a , a' , b , b' , c ; x , y ) = Sum_{m,n=0 to infinity} [ (a)_m (a')_n (b)_m (b')_n / (c)_m+n ] ( x^m / m! ) ( yⁿ /n! )

STACK	INPUTS	OUTPUTS
Level 7	a	/
Level 6	a'	/
Level 5	b	/
Level 4	b'	/
Level 3	c	/
Level 2	x	/
Level 1	y	F₃(a,b,b',c,c';x,y)

Example:

2 ENTER
3    ENTER
4    ENTER
6 ENTER
7 ENTER
   0.1 ENTER
0.2 APL3 -> F₃( 2 , 3 , 4 , 6 , 7 , 0.1 , 0.2 ) = 1.97231096375

Note:

-The series is convergent if Max { | x | , | y | } < 1

d) Appell Function F₄

-The Appell hypergeometric function F₄ is defined by:

F₄( a , b , c , c' ; x , y ) = Sum_{m,n=0
to infinity} [ (a)_m+n (b)_m+n / ( (c)_m (c')_n) ] ( x^m / m! ) ( yⁿ /n! )

STACK	INPUTS	OUTPUTS
Level 6	a	/
Level 5	b	/
Level 4	c	/
Level 3	c'	/
Level 2	x	/
Level 1	y	F₁(a,b,b',c;x,y)

Example:

2 ENTER
3    ENTER
4    ENTER
7 ENTER
   0.1 ENTER
0.2 APL4 -> F₄( 2 , 3 , 4 , 7 , 0.1 , 0.2 ) = 1.45948287121

Notes:

-The series is convergent if | x |^1/2 + | y |^1/2 < 1

-All these double series are special cases of "Kampé de Fériet" functions and may be computed by KdF

11°) Kampé de Fériet Functions

-In KAMP sub-directory, we have KdF functions

-Kampé de Fériet function is the generalized hypergeometric function of 2 variables.
-Given 6 non-negative integers: A , B , C , P , Q , S , it is defined by the double series:

           _{A B
C}   /    a₁ ,........, a_A ; b₁ ,........, b_B ; c₁ ,........, c_C ;              \
        F_{P Q}_S   |                                                                             x , y    |      =    SUM_m=0,1,2,... SUM_n=0,1,2,..    K_m,n x^m yⁿ / ( m! n! )
                      \    p₁ ,........, p_P ; q₁ ,........, q_Q ; s₁ ,........, s_S ;               /

where K_m,n = [ (a₁)_m+n ......... (a_A)_m+n (b₁)_m ......... (b_B)_m (c₁)_n ......... (c_C)_n ] / [ (p₁)_m+n ......... (p_P)_m+n (q₁)_m ......... (q_Q)_m (s₁)_n ......... (s_S)_n ]

and (a)_n = Pochhammer's symbol: (a)₀ = 1
(a)_n = a(a+1)(a+2) ...... (a+n-1) if n > 0

STACK	INPUTS	OUTPUTS
Level 2	x	/
Level 1	y	KdF(x,y)

Example:

A = 2 B = 3 C = 2 a₁ = 1.2 a₂ = 2.3 , b₁ = 1.1 b₂ = 0.7 b₃ = 1.4 , c₁ = 1.6 c₂ = 1.5

P = 3 Q = 2 S = 4 p₁ = 2 p₂ = 3 p₃ = 4 , q₁ = 5 q₂ = 6 , s₁ = 7 s₂ = 8 s₃ = 9 s₄ = 10.1

-Calculate F(5,7)

1°) Store the 6 constants A , B , C , P , Q , S in the corresponding variables
2°) Store the 16 other constants as a vector in'V' variables
3°)
5 ENTER
7 KdF -> F(5;7) = 1.02246670867

Note:

-This program doesn't check if the series is convergent ( press ON to stop the loops if they are infinte... or too slow ! )

12°) Lauricella Functions

a) Lauricella Function F_A

-This function is defined by

F_A( a ; b₁ , b₂ , b₃ ; c₁ , c₂ , c₃ ; x , y , z ) = SUM_{m,n,p=0
to infinity} (a)_m+n+p (b₁)_m (b₂)_n (b₃)_p / ((c₁)_m (c₂)_n (c₃)_p) x^myⁿz^p / (m! n! p!)

where (a)_k = a(a+1)(a+2) ....... (a+k-1) if k > 0 and (a)₀ = 1 ( Pochhammer symbol )

STACK	INPUTS	OUTPUTS
Level 10	a	/
Level 9	b₁	/
Level 8	b₂	/
Level 7	b₃	/
Level 6	c₁	/
Level 5	c₂	/
Level 4	c₃	/
Level 3	x	/
Level 2	y	/
Level 1	z	F_A(x,y,z)

Example: With a = 0.7 , b₁ = 0.1 b₂ = 0.2 b₃ = 0.3 , c₁ = 2.4 c₂ = 2.5 c₃ = 2.6 ,

F_A( 0.11 , 0.12 , 0.13 ) = ?

0.7 ENTER
0.1 ENTER
0.2 ENTER
0.3 ENTER
2.4 ENTER
2.5 ENTER
2.6 ENTER

0.11 ENTER
0.12 ENTER
0.13 LCFA -> F_A( 0.11 , 0.12 , 0.13 ) = 1.02157159441

Notes:

-With an HP48GX, execution time = 16 seconds
-The series is convergent if | x | + | y | + | z | < 1

b) Lauricella Function F_B

-This function is defined by

F_B( a₁ , a₂ , a₃ ; b₁ , b₂ , b₃ ; c ; x , y , z ) = SUM_{m,n,p=0 to infinity} (a₁)_m (a₂)_n (a₃)_p (b₁)_m (b₂)_n (b₃)_p / (c)_m+n+p x^myⁿz^p / (m! n! p!)

STACK	INPUTS	OUTPUTS
Level 10	a₁	/
Level 9	a₂	/
Level 8	a₃	/
Level 7	b₁	/
Level 6	b₂	/
Level 5	b₃	/
Level 4	c	/
Level 3	x	/
Level 2	y	/
Level 1	z	F_B(x,y,z)

Example: With a₁ = 0.1 a₂ = 0.2 a₃ = 0.3 , b₁ = 0.4 b₂ = 0.5 b₃ = 0.6 , c = 2.7

F_B( 0.41 , 0.42 , 0.43 ) = ?

0.1 ENTER
0.2 ENTER
0.3 ENTER
0.4 ENTER
0.5 ENTER
0.6 ENTER
2.7 ENTER

0.41 ENTER
0.42 ENTER
0.43 LCFB -> F_B( 0.41 , 0.42 , 0.43 ) = 1.05774803114

Note:

-The series is convergent if Max { | x | , | y | , | z | } < 1

c) Lauricella Function F_C

-This function is defined by

F_C( a ; b ; c₁ , c₂ , c₃ ; x , y , z ) = SUM_{m,n,p=0 to infinity} (a)_m+n+p (b)_m+n+p / ( (c₁)_m(c₂)_n(c₃)_p ) x^myⁿz^p / (m! n! p!)

STACK	INPUTS	OUTPUTS
Level 8	a	/
Level 7	b	/
Level 6	c₁	/
Level 5	c₂	/
Level 4	c₃	/
Level 3	x	/
Level 2	y	/
Level 1	z	F_C(x,y,z)

Example: With a = 0.1 , b = 0.2 , c₁ = 2.7 c₂ = 2.8 c₃ = 2.9

F_C( 0.10 , 0.11 , 0.12 ) = ?

0.1 ENTER
0.2 ENTER
2.7 ENTER
2.8 ENTER
2.9 ENTER

0.10 ENTER
0.11 ENTER
0.12 LCFC -> F_C( 0.10 , 0.11 , 0.12 ) = 1.00255428904

Note:

-The series is convergent if | x |^1/2 + | y |^1/2 + | z |^1/2 < 1

d) Lauricella Function F_D

-This function is defined by

F_D( a ; b₁ , b₂ , b₃ ; c ; x , y , z ) = SUM_{m,n,p=0 to infinity} (a)_m+n+p (b₁)_m (b₂)_n (b₃)_p / (c)_m+n+p x^myⁿz^p / (m! n! p!)

STACK	INPUTS	OUTPUTS
Level 8	a	/
Level 7	b₁	/
Level 6	b₂	/
Level 5	b₃	/
Level 4	c	/
Level 3	x	/
Level 2	y	/
Level 1	z	F_D(x,y,z)

Example: With a = 0.1 , b₁ = 0.2 b₂ = 0.3 b₃ = 0.4 , c = 2.7

F_D( 0.21 , 0.31 , 0.41 ) = ?

0.1 ENTER
0.2 ENTER
0.3 ENTER
0.4 ENTER
2.7 ENTER

0.21 ENTER
0.31 ENTER
0.41 LCFD -> F_D( 0.21 , 0.31 , 0.41 ) = 1.01234537728

Notes:

-The series is convergent if Max { | x | , | y | , | z | } < 1
-All these triple series are special cases of a generalized hypergeometric series of 3 variables: Srivastava functions

13°) Srivastava Functions

-Srivastava's Functions generalize the hypergeometric functions of 3 variables x , y , z.

-Given 14 non-negative integers A B C D E F G H I J K L M N and

A numbers a₁ ..... a_A B numbers b₁ ..... b_B C numbers c₁ ... c_C D numbers d₁ .... d_D E numbers e₁ .... e_E F numbers f₁ .... f_F G numbers g₁....g_G
H numbers h₁ ..... h_H I numbers i₁ ..... i_I J numbers j₁ ... j_J K numbers k₁ .... k_K L numbers l₁ .... l_L M numbers m₁ ... m_M N numbers n₁ ... n_N

the following program computes

F ^A_H^B_I^C_J^D_K^E_L^F_M^G_N(^{a1 ... aA ;}_{h1 ... hH ;}^{b1 ...
bB ;}_{i1 ... iI ;}^{c1 ... cC ;}_{j1 ...
jJ ;}^{d1 ... dD ;}_{k1 ... kK ;}^{e1 ...
eE ;}_{l1 ...lL ;}^{f1 ... fF ;}_{m1 ... mM
;}^{g1... gG ;}_{n1 ... nN} ; x , y , z )

= SUM_{p,q,r=0 to infinity} u_p,q,r x^p y^q z^r / ( p! q! r! )

where u_p,q,r = [ (a₁)_p+q+r ...... (a_A)_p+q+r (b₁)_p+q ...... (b_B)_p+q (c₁)_q+r ...... (c_C)_q+r (d₁)_p+r ...... (d_D)_p+r (e₁)_p ...... (e_E)_p (f₁)_q ...... (f_F)_q (g₁)_r ...... (g_G)_r ] /
[ (h₁)_p+q+r ...... (h_H)_p+q+r (i₁)_p+q ...... (i_I)_p+q (j₁)_q+r ...... (j_J)_q+r (k₁)_p+r ...... (k_K)_p+r (l₁)_p ...... (l_L)_p (m₁)_q ...... (m_M)_q (n₁)_r ...... (n_N)_r ]

and (t)_n = Pochhammer's symbol:

(t)₀ = 1
(t)_n = t(t+1)(t+2) ...... (t+n-1) if n > 0

STACK	INPUTS	OUTPUTS
Level 3	x	/
Level 2	y	/
Level 1	z	F(x,y,z)

Example: With A = B = ........................... = M = 1 & N = 2

a₁ = sqrt(2) b₁ = sqrt(3) c₁ = sqrt(5) d₁ = sqrt(6) e₁ = sqrt(7) f₁ = sqrt(8) g₁ = sqrt(10)
h₁ = sqrt(11) i₁ = sqrt(12) j₁ = sqrt(13) k₁ = sqrt(14) l₁ = sqrt(15) m₁ = sqrt(17) n₁ = sqrt(18) n₂ = sqrt(19)

and x = 0.1 y = 0.2 z = 0.3

1°) Store 1 in all the 14 variables 'A' 'B' ....... 'N'
2°) Place the 15 constants sqrt(2) sqrt(3) ..... sqrt(19) in the stack & 15 ->ARRY and store this vector in 'V' variable
3°)

0.1 ENTER
0.2 ENTER
0.3 SHRI -> F( 0.1 , 0.2 , 0.3 ) = 1.03769933592

Notes:

-With an HP48GX, execution time = 37 seconds in the example above
-This program does not check if the series is convergent or not.

14°) Associated Legendre Functions

ALF1 ALF3 ALF5 compute the associated Legendre functions of the 1st kind P_n^m(x)

ALF1 if | x | < 1
ALF3 if | 1-x | < 2
ALF5 if | x | > 1

-If | x | < 1 , they use what is called "Associated Legendre function of the first kind of type 2" in reference [2]
-If | x | > 1 , ---------------------- "Associated Legendre function of the first kind of type 3 ---------------
-Cf http://functions.wolfram.com/ to see all the formullae...

-Likewise, ALF2 ALF4 ALF6 compute the associated Legendre functions of the 2nd kind Q_n^m(x)

ALF2 if | x | < 1
ALF4 if | 1-x | < 2
ALF6 if | x | > 1

STACK	INPUTS	OUTPUTS
Level 3	m	/
Level 2	n	/
Level 1	x	P_n^m(x ) or Q_n^m(x)

Examples:

1.2 ENTER
2.3 ENTER
0.7 ALF1 -> P_2.3^1.2(0.7) = -1.87245992968

1.2 ENTER
2.3 ENTER
0.7 ALF2 -> Q_2.3^1.2(0.7) = 1.17884238034

1.2 ENTER
2.3 ENTER
1.4 ALF3 -> P_2.3^1.2(1.4) = 6.13299598995

1.2 ENTER
2.3 ENTER
1.4 ALF4 -> Q_2.3^1.2(1.4) = -0.238771340462 - 0.173477533315 i ( complex number )

1.2 ENTER
2.3 ENTER
2.4 ALF5 -> P_2.3^1.2(2.4) = 29.6468224384

1.2 ENTER
2.3 ENTER
2.4 ALF6 -> Q_2.3^1.2(1.4) = -0.0243521404172 - 0.0176928656611 i ( complex number )

1 ENTER
2 ENTER
3 ALF6 -> Q₂¹(3) = -0.016511473615

-If the m & n-values are integers, with 0 <= m <= n we can use the following formulae, employed by PMN:

(n-m) P_n^m(x) = x (2n-1) P_n-1^m(x) - (n+m-1) P_n-2^m(x)

P_m^m(x) = (-1)^m (2m-1)!! ( 1-x² )^m/2 if | x | < 1 where (2m-1)!! = (2m-1)(2m-3)(2m-5).......5.3.1
P_m^m(x) = (2m-1)!! ( x²-1 )^m/2 if | x | > 1

-If m = 0 , we have Legendre polynomials

Example:

4 ENTER
   7 ENTER
   3 PMN ->    P₇⁴(3) = 37920960

-We also have the relations used by PQMN

P_n^m(x) = ( x²-1 )^m/2 d^mP_n(x)/dx^m if | x | > 1 P_n^m(x) = (-1)^m ( 1- x² )^m/2 d^mP_n(x)/dx^m if | x | < 1
Q_n^m(x) = ( x²-1 )^m/2 d^mQ_n(x)/dx^m if | x | > 1 Q_n^m(x) = (-1)^m ( 1- x² )^m/2 d^mQ_n(x)/dx^m if | x | < 1

-if m = 0 P_n^m(x) = P_n(x) = Legendre polynomials
and Q_n^m(x) = Q_n(x) with Q₀(x) = 0.5 Ln | (1+x)/(1-x) | ; Q₁(x) = x/2 Ln | (1+x)/(1-x) | - 1 and n.Q_n(x) = (2n-1).x.Q_n-1(x) - (n-1).Q_n-2(x)

-We have employed the recurrence relations: P_n^m+1(x) = | x²-1 |^-1/2 [ (n-m).x. P_n^m(x) - (n+m). P_n-1^m(x) ] ( the same relation holds for Q_n^m(x) )

Important note:

-If | x | is significantly greater than 1 , Q_n^m(x) is obtained very inaccurately ( the recurrence relation is unstable in this case )

STACK	INPUTS	OUTPUTS
Level 3	m	/
Level 2	n	P_n^m(x )
Level 1	x	Q_n^m(x)

Examples:

    4      ENTER
    7      ENTER    P₇⁴(0.6) = 715.309056
0.6     PQMN yields   Q₇⁴(0.6) = -1011.1718046

-Similarly: P₇⁴(1.2) = 6327.21196829 & Q₇⁴(1.2) = 82.1210783385

-With x = 1.2 Q₇⁴(x) is still acceptable ( ALF6 produces 82.1210783254 )
but for x = 3 , the errors are too large and ALF6 is much better

15°) Bessel Functions

a) Bessel & modified Bessel Functions of the 1st & 2nd kind

JNX computes Bessel function of the 1st kind J_n(x) = (x/2)ⁿ [ 1/Gam(n+1) + (-x²/4)¹/ (1! Gam(n+2) ) + .... + (-x²/4)^k/ (k! Gam(n+k+1) ) + .... ]

STACK	INPUTS	OUTPUTS
Level 2	n	/
Level 1	x	J_n(x)

Example:

0.7   ENTER
   1.9   JNX   yields   J_0.7(1.9) = 0.584978103022

YNX returns Bessel function of the 2nd kind:    Y_n(x) = ( J_n(x) cos(n(pi)) - J_-n(x) ) / sin(n(pi))

STACK	INPUTS	OUTPUTS
Level 2	n	/
Level 1	x	Y_n(x)

Example:

1.2 ENTER
1.4 YNX -> Y_1.2(1.4) = -0.593492008383

Note:

YNX1 may also be used instead of YNX
It calls UABX as a subroutine

INX gives the modified Bessel function of the 1st kind I_n(x) = (x/2)ⁿ [ 1/Gam(n+1) + (x²/4)¹/ (1! Gam(n+2) ) + .... + (x²/4)^k/ (k! Gam(n+k+1) ) + .... ]

STACK	INPUTS	OUTPUTS
Level 2	n	/
Level 1	x	I_n(x)

Example:

0.7 ENTER
1.9 INX yields I_0.7(1.9) = 1.72763060315

KNX calculates the modified Bessel function of the 2nd kind K_n(x) = (pi/2) ( I_-n(x) - I_n(x) ) / sin(n(pi))

STACK	INPUTS	OUTPUTS
Level 2	n	/
Level 1	x	K_n(x)

Example:

1.2 ENTER
1.4 KNX -> K_1.2(1.4) = 0.361241917679

Note:

KNX1 may also be used instead of KNX
It calls UABX as a subroutine

-The programs above also work with complex arguments.
-They employ the hypergeometric functions
-So, the precision is not very good for large x-values ( it may even become meaningless results )

-The following program JYNX employs a continued fraction
-They also work well with large x-values, but the arguments must remain real

Formulae:

With p + i.q = -1/(2x) + i + (i/x) [ ( 0.5² - n² )/( 2x + 2i + ( 1.5² - n² )/( 2x + 4i + ..... ) ) ]
and g_n = -1/(((2n + 2)/x) - 1/(((2n + 4)/x) - ..... ))

J_n(x) = sign(D) [ ( 2q/(x.Pi) ) / ( q² + ( p - g_n - n/x )² ) ]^1/2 where D = the denominator of the second continued fraction
Y_n(x) = [ ( p - g_n - n/x )/q ] J_n(x)

STACK	INPUTS	OUTPUTS
Level 2	n	J_n(x)
Level 1	x	Y_n(x)

Example:

3.14 ENTER J_3.14(100) = 0.079535723253
100 JYNX -> Y_3.14(100) = 0.00658232732688653

ITJ computes the integral §₀^x J_n(t).dt = 2 ( J_n+1(x) + J_n+3(x) + J_n+5(x) + ........ )

STACK	INPUTS	OUTPUTS
Level 2	n	/
Level 1	x	§₀^x J_n(t).dt

Example:

1 ENTER
     3 ITJ ->       §₀³ J₁(x).dx = 1.2600519549

     b) Spherical Bessel Functions

-Spherical Bessel functions of the first kind are defined by   j_n(x) = (Pi/(2x))^1/2 J_n+1/2(x)   where n is an integer
-We can also compute them by the following formula:

           j_n-1(x) = j_n(x) (2n+1)/x - j_n+1(x)   this recurrence relation is used , starting with j_m(x) = 1 , j_m+1(x) = 0 for some large m

-Then, the results are normalized by the sum     1 j₀²(x) + 3 j₁²(x) + 5 j₂²(x) + ....... = 1

STACK	INPUTS	OUTPUTS
Level 2	n > 0	/
Level 1	x	j_n(x)

Example:

2 PI SB1 yields j₂(Pi) = 0.303963550927

-Spherical Bessel functions of the second kind are defined by y_n(x) = (Pi/(2x))^1/2 Y_n+1/2(x) where n is an integer
-We also have::

y_n+1(x) = y_n(x) (2n+1)/x - y_n-1(x)

y₀(x) = -(cos x)/x , y₁(x) = -(cos x)/x² - (sin x)/x

STACK	INPUTS	OUTPUTS
Level 2	n > 0	/
Level 1	x	y_n(x)

Example:

2 PI SB2 -> -0.221555282885

c) Kelvin Functions

-Kelvin functions of the 1st kind are defined by:

ber_n(x) = (x/2)ⁿ Sum_{k=0,1,2,......} cos ((3n/4+k/2)180°) (x²/4)^k / ( k! Gam(n+k+1) ) where Gam = Euler's Gamma function

bei_n(x) = (x/2)ⁿ Sum_{k=0,1,2,......} sin ((3n/4+k/2)180°) (x²/4)^k / ( k! Gam(n+k+1) )

STACK	INPUTS	OUTPUTS
Level 2	n	ber_n(x)
Level 1	x	bei_n(x)

Example:

   2   SQRT    ber_sqrt(2)(PI) = -0.674095952864
PI   KLV1 ->      bei_sqrt(2)(PI) = -1.59735721106

-Kelvin functions of the 2nd kind are defined by:

ker_n(x) = - (PI/2) [ bei_n(x) - (ber_-n(x)/sin(n.180°)) + (ber_n(x)/tan(n.180°)) ] ( cf reference [2] for other formulae if n is an integer )

kei_n(x) = (PI/2) [ ber_n(x) + (bei_-n(x)/sin(n.180°)) - (bei_n(x)/tan(n.180°)) ]

STACK	INPUTS	OUTPUTS
Level 2	n	ker_n(x)
Level 1	x	kei_n(x)

Example:

2 SQRT ker_sqrt(2)(PI) = 0.0259018942636
PI KLV2 -> kei_sqrt(2)(PI) = 0.089242864903

Note:

S1 & S2 are subroutines called by the Bessel function programs

16°) Orthogonal Polynomials & Functions

a) Laguerre Polynomial

Formulae: L_n(x) = (2n-1-x).L_n-1(x) - (n-1)².L_n-2(x) ; L₀(x) = 1 ; L₁(x) = 1 - x

STACK	INPUTS	OUTPUTS
Level 2	n >= 0	/
Level 1	x	L_n(x)

Example:

7 ENTER
3.14 LAG -> L₇ (3.14) = -4932.43994869

Note:

-Some authors divide L_n (x) by n!

b) Generalized Laguerre Polynomials & Functions

Formulae: L₀^(a) (x) = 1 ; L₁^(a) (x) = a+1-x ; n.L_n^(a) (x) = (2.n+a-1-x).L_n-1^(a) (x) - (n+a-1).L_n-2^(a) (x)

STACK	INPUTS	OUTPUTS
Level 3	a	/
Level 2	n >= 0	/
Level 1	x	L_n^a(x)

Example:

1.4 ENTER
7 ENTER
PI PANX -> L₇^(1.4)(Pi) = 1.68889351366

-For a non integer n-value, we can use L_n^(a)(x) = [ Gam(n+a+1) / Gam(n+1) / Gam(a+1) ] M(-n,a+1,x)

where Gam = Gamma function
and M = Kummer's function

STACK	INPUTS	OUTPUTS
Level 3	a	/
Level 2	n	/
Level 1	x	L_n^a(x)

Example:

2   SQRT
7   SQRT
PI 1/X LANX -> L_sqrt(7)^(sqrt(2))(1/PI) = 3.62560918475

     c) Legendre Polynomials

Formulae:      n.P_n(x) = (2n-1).x.P_n-1(x) - (n-1).P_n-2(x) ; P₀(x) = 1 ; P₁(x) = x

STACK	INPUTS	OUTPUTS
Level 2	n >= 0	/
Level 1	x	P_n(x)

Example:

7 ENTER
3.14 LEG -> P₇(3.14) = 68076.394314

d) Jacobi Polynomials & Functions

Formulae: P₀^(a;b) (x) = 1 ; P₁^(a;b) (x) = (a-b)/2 + x.(a+b+2)/2

2n(n+a+b)(2n+a+b-2).P_n^(a;b) (x) = [ (2n+a+b-1).(a²-b²) + x.(2n+a+b-2)(2n+a+b-1)(2n+a+b) ] P_n-1^(a;b) (x)
- 2(n+a-1)(n+b-1)(2n+a+b) P_n-2^(a;b) (x)

STACK	INPUTS	OUTPUTS
Level 4	a	/
Level 3	b	/
Level 2	n >= 0	/
Level 1	x	P_n^(a;b) (x)

Example:

1.4 ENTER
1.7 ENTER
7 ENTER
PI 1/X JCP -> P₇^(1.4;1.7)(1/Pi) = -0.32223442011

Note:

-The hypergeometric function also allows to evaluate Jacobi's functions even when n is not an integer:

P_n^(a;b) (x) = [ Gam(a+n+1) / Gam(a+1) / Gam(n+1) ] F ( -n , a+b+n+1 , a+1 , (1-x)/2 )

STACK	INPUTS	OUTPUTS
Level 4	a	/
Level 3	b	/
Level 2	n	/
Level 1	x	P_n^(a;b) (x)

Example:

   2 SQRT
   3 SQRT
   7 SQRT
PI 1/X JCF -> P_sqrt(7)^{(sqrt(2),sqrt(3))}(1/PI) = -0.887945297368

     e) Hermite Polynomials & Functions

-Hermite polynomials are defined by:

      H_n(x) = 2x.H_n-1(x) - 2(n-1).H_n-2(x) ; H₀(x) = 1 ; H₁(x) = 2x

STACK	INPUTS	OUTPUTS
Level 2	n	/
Level 1	x	P_n(x)

Example:

7 ENTER
3.14 HMT -> H₇ (3.14) = 73726.2432426

Note:

-With non-integer n-values, we have H_n(x) = 2ⁿ sqrt(PI) [ (1/Gam((1-n)/2)) M(-n/2,1/2,x²) - ( 2.x / Gam(-n/2) ) M((1-n)/2,3/2,x²) ]

where Gam = Gamma function
and M = Kummer's function

STACK	INPUTS	OUTPUTS
Level 2	n	/
Level 1	x	P_n(x)

Example:

   0.4   ENTER
    PI   1/X HMTF ->    H_0.4(1/PI) = 1.01032199692

     f) Chebyshev Polynomials & Functions

-Chebyshev polynomials are defined by:

Formulae:    T_n(x) = 2x.T_n-1(x) - T_n-2(x) ; T₀(x) = 1 ; T₁(x) = x    ( first kind )

U_n(x) = 2x.U_n-1(x) - U_n-2(x) ; U₀(x) = 1 ; U₁(x) = 2x    ( second kind )

STACK	INPUTS	OUTPUTS
Level 2	n >= 0	T_n(x)
Level 1	x	U_n(x)

Examples:

7 ENTER T₇ (0.314) = -0.786900700394
0.314 CHB -> U₇ (0.314) = -0.582815680411

-If n is not an integer, we can also use: T_n(x) = Cos ( n Cos^-1 x ) & U_n(x) = Sin ( (n+1) Cos^-1 x ) / sqrt ( 1-x² )

STACK	INPUTS	OUTPUTS
Level 2	n	/
Level 1	x	T_n(x) or U_n(x)

Examples:

   7.12    ENTER
0.314 TNX ->   T_7.12 (0.314) = -0.870362214734

   7.12    ENTER
0.314 UNX ->   U_7.12 (0.314) = -0.707508162579

     g) Ultraspherical Polynomials & Functions

-The polynomials are defined by:

   C₀^(a) (x) = 1 ;   C₁^(a) (x) = 2.a.x   ;    (n+1).C_n+1^(a) (x) = 2.(n+a).x.C_n^(a) (x) - (n+2a-1).C_n-1^(a) (x)      if a # 0      C_n⁽⁰⁾ (x) = (2/n).T_n(x)

STACK	INPUTS	OUTPUTS
Level 3	a	/
Level 2	n >= 0	/
Level 1	x	C_n^(a)(x)

Example:

1.5 ENTER^
7 PI 1/X USP -> C₇^(1.5)(1/Pi) = -0.989046385317

Notes:

-If n is not an integer, Ultraspherical functions are defined as:

C_n⁽⁰⁾(x) = (2/n) T_n(x) where T_n(x) = Chebyshev function, first kind and if a#0:
C_n^(a)(x) = [ Gam(n+2a) / Gam(n+1) / Gam(2a) ] F( -n , n+2a , n+1/2 , (1-x)/2 ) where F is the hypergeometric function and Gam = Gamma function

STACK	INPUTS	OUTPUTS
Level 3	a	/
Level 2	n	/
Level 1	x	C_n^(a)(x)

Example:

2 SQRT
7 SQRT
PI 1/X USF -> C_sqrt(7)^(sqrt(2))(1/PI) = -1.57546640046

17°) Elliptic Functions

a) WeierStrass

-The Weierstrass Elliptic Function P(x;g₂;g₃) satisfies the differential equation: (P')² = 4.P³ - g₂.P - g₃
-WEF calculates P(x;g₂;g₃) by a Laurent series:

P(x;g₂;g₃) = x^-2 + c₂.x² + c₃.x⁴ + ...... + c_n.x^2n-2 + ....

where c₂ = g₂/20 ; c₃ = g₃/28 and c_n = 3 ( c₂. c_n-2 + c₃. c_n-3 + ....... + c_n-2. c₂ ) / (( 2n+1 )( n-3 )) ( n > 3 )

STACK	INPUT	OUTPUT
Level 1	x	P(x;g₂;g₃)

Example: Calculate P(x;g₂;g₃) for x = 0.6 g₂ = 0.9 g₃ = 1.4

-First, store g₂ & g₃ in the corresponding variables: 0.9 'G2' STO 1.4 'G3' STO

0.6 WEF -> P( 0.6 ; 0.9 ; 1.4 ) = 2.8005007841

Duplication Formula:

P(2x) = -2 P(x) + ( 6 P²(x) - g₂/2 )² / ( 4 ( 4 P³(x) - g₂ P(x) - g₃ ) )

-This formula is used recursively ( n times ) to obtain P(2ⁿx) if we know P(x) ( n = 1 ; 2 ; 3 ; ..... )

STACK	INPUTS	OUTPUTS
Level 2	P(x)	/
Level 1	n > 0	P(2ⁿx)

Example: Calculate P(x;g₂;g₃) for x = 4.8 g₂ = 0.9 g₃ = 1.4

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

( without changing 'G2' & 'G3' )

2.8005007841 ENTER
3 DUPW -> P(4.8;0.9;1.4) = 1.9547041716

b) Jacobian Elliptic Functions

JEF calculates sn ( u | m ) , cn ( u | m ) & dn ( u | m )

STACK	INPUTS	OUTPUTS
Level 3	/	sn ( u \| m )
Level 2	u	cn ( u \| m )
Level 1	m	dn ( u \| m )

Examples:

1- Evaluate sn ( 0.7 | 0.3 ) cn ( 0.7 | 0.3 ) dn ( 0.7 | 0.3 )

sn ( 0.7 | 0.3 ) = 0.63230477631
0.3 ENTER cn ( 0.7 | 0.3 ) = 0.774719736327
0.7 JEF -> dn ( 0.7 | 0.3 ) = 0.938113639679

2-Likewise,

                sn ( 0.7 | 1 ) = 0.604367777117          sn ( 0.7 | 2 ) = 0.564297007561        sn ( 0.7 | -3 ) = 0.759113420483
                cn ( 0.7 | 1 ) = 0.796705459993          cn ( 0.7 | 2 ) = 0.825571854710        cn ( 0.7 | -3 ) = 0.650958381789
                dn ( 0.7 | 1 ) = 0.796705459993          dn ( 0.7 | 2 ) = 0.602609139091        dn ( 0.7 | -3 ) =1.65189574594

Note:

JEF also works with complex numbers.

c) Carlson Elliptic Integrals

  RC is defined by

R_C(x;y) = R_F(x;y;y)    if y > 0
R_C(x;y) = (x/(x-y))^1/2 R_C(x-y;-y)      if y < 0

STACK	INPUTS	OUTPUTS
Level 2	x	/
Level 1	y	R_C(x;y)

Examples:

1 ENTER
3 RC -> R_C(1;3) = 0.675510858854

1 ENTER
-3 RC -> R_C(1;-3) = 0.274653072166

Notes:

RC also works with complex arguments
It's actually a degenerate form of RF:

-Carlson has given a new definition of a standard elliptic integral of the first kind:

R_F(x;y;z) = (1/2) §₀^infinity ( ( t + x ).( t + y ).( t + z ) )^-1/2 dt with x , y , z non-negative and at most one is zero

RF uses the following property:

R_F(x;y;z) = R_F((x+L)/4;(y+L)/4;(z+L)/4) where L = x^1/2y^1/2 + x^1/2z^1/2 + y^1/2z^1/2

-This transformation is performed until x , y , z are close enough to apply R_F(x;y;z) = µ^-1/2 with µ = (x+y+z)/3 ( we have R_F(x;x;x) = x^-1/2 )

STACK	INPUTS	OUTPUTS
Level 3	x	/
Level 2	y	/
Level 1	z	R_F(x;y;z)

Example:

    2 ENTER
    3 ENTER
    4 RF ->    R_F(2;3;4) = 0.584082841671

Note:

RC also works with complex arguments

-The elliptic integral of the third kind RJ is defined by

R_J(x;y;z;p) = (3/2) §₀^infinity ( t + p )^-1 ( ( t + x ).( t + y ).( t + z ) )^-1/2 dt with x , y , z non-negative and at most one is zero p > 0

-We have R_J(x,x,x,x) = x^-3/2

RJ applies R_J(x;y;z;p) = 3 Sum_{n=0,1,2,....,k} R_F(a_n,b_n,b_n)/4ⁿ + 1/(4^k+1µ^3/2)

where x₀ = x , y₀ = y , z₀ = z , p₀ = p ; x_n+1 = ( x_n+ L_n )/4 , y_n+1 = ( y_n+ L_n )/4 , z_n+1 = ( z_n + L_n )/4 , p_n+1 = ( p_n +L_n )/4

with L_n = x_n^1/2y_n^1/2 + x_n^1/2z_n^1/2 + y_n^1/2z_n^1/2
a_n = ( p_n( x_n^1/2 + y_n^1/2 + z_n^1/2 ) + ( x_ny_nz_n )^1/2 )² ; b_n = p_n ( p_n + L_n )²

and µ = (x_k+1+y_k+1+z_k+1+2p_k+1)/5

-The iterations stop when x_n , y_n , z_n , p_n are close enough to produce a good accuracy.
-This program also computes the Cauchy principal value of the integral if p < 0

STACK	INPUTS	OUTPUTS
Level 4	x	/
Level 3	y	/
Level 2	z	/
Level 1	p # 0	R_J(x;y;z;p)

Examples:

     1 ENTER
     2 ENTER
     3 ENTER
     4 RJ ->    R_J(1;2;3;4) = 0.239848099749

    1 ENTER
     2 ENTER
     3 ENTER
    -4    RJ        >>>> R_J(1;2;3;-4) = -0.23786769473

Notes:

RJ also work if x , y , z are complex numbers if p > 0
But if p < 0 , x , y , z must be real numbers

RFZ computes R_F( a ; b+i.c ; b-i.c ) where a , b , c are real numbers

STACK	INPUTS	OUTPUTS
Level 3	a	/
Level 2	b	/
Level 1	c	R_F(a;b+i.c;b-i.c)

Example:

      2 ENTER
      3 ENTER
      4 RFZ -> R_F ( 2 , 3+4i , 3-4i ) = 0.543421944621

-This result could be obtained by RF but more slowly
-Likewise, RJZ calculates R_J( a ; b+i.c ; b-i.c ; p ) where a , b , c , p are real numbers ( p # 0 )

STACK	INPUTS	OUTPUTS
Level 4	a	/
Level 3	b	/
Level 2	c	/
Level 1	p # 0	R_J(a;b+i.c;b-i.c;p)

Example:

    1 ENTER
    2 ENTER
    3 ENTER
    4 RJZ -> R_J ( 1 , 2+3i , 2-3i , 4 ) = 0.20564414055

Note:

-Unfortunately, RJZ doesn't always work if p < 0
-A formula is used to change p into a positive number but it sometimes fails:

R_J ( 2 , 3+4i , 3-4i , -7 ) = -0.0939148073957 is correctly computed

but with R_J ( 2 , 3+4i , 3-4i , -6 ) , -6 becomes -0.125 so we cannot use RJZ ( the HP48 displays "Bad Argument Value" )

-Carlson has also defined a symmetric Elliptic Integral of the second kind:

R_G(x;y;z) = (1/4) §₀^infinity ( ( t + x ).( t + y ).( t + z ) )^-1/2 .( x/(t+x) + y/(t+y) + z/(t+z) ) t.dt

And we have: 2.R_G(x;y;z) = z.R_F(x;y;z) - (x-z)(y-z)/3 R_D(x;y;z) + ( x.y/z )^1/2

STACK	INPUTS	OUTPUTS
Level 3	x	/
Level 2	y	/
Level 1	z	R_G(x;y;z)

Example:

    2 ENTER
    3 ENTER
    4 RG ->    R_G(2;3;4) = 1.72550302806

     d) Surface Area of an Ellipsoid

-The symmetric elliptic integral of the second kind may also be used to compute the area of an ellipsoid:

Equation of an ellipsoid: x²/a² + y²/b² +x²/c² = 1

Area = 4.PI.R_G( a²b² , a²c² , b²c² )

STACK	INPUTS	OUTPUTS
Level 3	a	/
Level 2	b	/
Level 1	c	Area

Example:

    2 ENTER
    4 ENTER
    9 SAE -> Area = 283.427384268

Note:

-Variable 'eps' contains a small number to stop the iterations for computing elliptic integrals ( eps = 10^-11 )

18°) Miscellaneous Functions

     a) Airy Functions

-They are defined by

    Ai(x) =   f(x) - g(x)                             with            f(x) = [ 3^-2/3 / Gamma(2/3) ] [ 1 + x³/3! + ( 1 . 4 x⁶ )/6! + ( 1 . 4 . 7 x⁹ )/9! + ..... ]
    Bi(x) = [ f(x) + g(x) ] sqrt(3)               and            g(x) = [ 3^-1/3 / Gamma(1/3) ] [ x + 2 x⁴/4! + ( 2 . 5 x⁷ )/7! + ( 2 . 5 . 8 x¹⁰ )/10! + ..... ]

-We can also use Hypergeometric functions

STACK	INPUTS	OUTPUTS
Level 2	/	Ai(x)
Level 1	x	Bi(x)

Example:

Ai(0.4) = 0.254742354289
0.4 AIBI ->   Bi(0.4) = 0.801773000017


     b) Anger & Weber Functions

Formulae:

J_n(x) = + (x/2) sin ( 90°n ) ₁F₂( 1 ; (3-n)/2 , (3+n)/2 ; -x²/4 ) / Gam((3-n)/2) / Gam ((3+n)/2) ( Anger )
+ cos ( 90°n ) ₁F₂( 1 ; (2-n)/2 , (2+n)/2 ; -x²/4 ) / Gam((2-n)/2) / Gam ((2+n)/2)

E_n(x) = - (x/2) cos ( 90°n ) ₁F₂( 1 ; (3-n)/2 , (3+n)/2 ; -x²/4 ) / Gam((3-n)/2) / Gam ((3+n)/2) ( Weber )
+ sin ( 90°n ) ₁F₂( 1 ; (2-n)/2 , (2+n)/2 ; -x²/4 ) / Gam((2-n)/2) / Gam ((2+n)/2)

STACK	INPUTS	OUTPUTS
Level 2	n	*J_n*(x)
Level 1	x	*E_n*(x)

Example:

   2   SQRT J_sqrt(2)(PI) = 0.366086558358
PI ANWEB ->    E_sqrt(2)(PI) = -0.315594384965


     c) Dawson Integral

-It is defined by F(x) = e^{-x^2} §₀^x e^{t^2} dt

DAW employs a continued fraction

STACK	INPUT	OUTPUT
Level 1	x	F(x)

Examples:

   1.94 DAW ->    F( 1.94 ) = 0.314057165471
100 DAW -> F(100) = 0.0050002500375


     d) Exponential Integrals E_n

-We have E_n(x) = §₁^+infinity t^-n e^-x.t dt ( x > 0 ; n = a positive integer )

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

-It may also be computed by a continued fraction

E_n(x) = exp (-x) ( 1/(x+n-1.n/(x+n+2-2(n+1)/(x+n+4- .... ))) )

-We also have: E₀(x) = (1/x).exp(-x) and E_n(0) = 1/(n-1) if n > 1

STACK	INPUTS	OUTPUTS
Level 2	n	/
Level 1	x	E_n(x)

Examples:

0 ENTER
1.4 ENX -> E₀(1.4) = 0.17614068853

2 ENTER
1.4 ENX -> E₂(1.4) = 0.0838899263411

    100 ENTER
     1.4 ENX -> E₁₀₀(1.4) = 0.00245580064836


     e) Exponential Integrals Ei

Ei(x) = §_-infinity^x e^t/t dt =     Ei(x) = C + ln(x) + Sum_n=1,2,..... xⁿ/(n.n!)    where C = 0.577215664902 is Euler's constant.

STACK	INPUTS	OUTPUTS
Level 1	x	E_i(x)

Example:

1.4 EI -> Ei(1.4) = 3.00720746416

Note:

Logarithmic Integral = Li(x) = Ei(Ln(x))    ( x > 1 )


     f) Sine & Cosine Integrals

Si(x) = §₀^x sin(t)/t dt      ;      Ci(x) = C + Ln(x) + §₀^x (cos(t)-1)/t dt

   Si(x) = Sum_{n=0,1,2,.....} (-1)ⁿ x²ⁿ⁺¹/((2n+1).(2n+1)!)
Ci(x) = C + ln(x) + Sum_n=1,2,..... (-1)ⁿ x²ⁿ/(2n.(2n)!)                 where C = 0.577215664902 is Euler's constant.

STACK	INPUTS	OUTPUTS
Level 1	x	Ci(x) or Si(x)

Example:

   1.4   SI     -> Si(1.4) = 1.25622673278
   1.4   CI    -> Ci(1.4) = 0.462006585095

Note:

-It also work with complex x-values
-But for large numbers, the precision becomes very small

-If x > 2 , it's preferable to use the following continued fraction:

   -Ci(x) + i.Si(x) - i.pi/2 = exp(-i.x) ( 1/(1+i.x-1²/(3+i.x-2²/(5+i.x-3²/(7+i.x- ..... )))) )

STACK	INPUTS	OUTPUTS
Level 3	/	Si(x) - Pi/2
Level 2	/	Ci(x)
Level 1	x	Si(x)

Examples:

Si(3)-pi/2 = 0.27785620121
Ci(3) = 0.119629786002
3 CISI -> Si(3) = 1.84865252801

Si(100)-pi/2 = -0.00857085990583
Ci(100) = -0.0051488251426
100 CISI -> Si(100) = 1.56222546689

Note: x must be a real number

     g) Hyperbolic Sine & Cosine Integrals

Definition:           Shi(x) = §₀^x sinh(t)/t dt      ;      Chi(x) = C + Ln(x) + §₀^x (cosh(t)-1)/t dt

Shi(x) = Sum_{n=0,1,2,.....} x²ⁿ⁺¹/((2n+1).(2n+1)!)
Chi(x)= C + ln(x) + Sum_n=1,2,..... x²ⁿ/(2n.(2n)!)    where C = 0.577215664902 is Euler's constant.

STACK	INPUTS	OUTPUTS
Level 1	x	Ci(x) or Si(x)

Example:

   1.4   SHI     -> Shi(1.4) = 1.56171338837
   1.4   CHI    -> Chi(1.4) = 1.44549407579

     h) Catalan Function

-Catalan numbers may be generalized to real ( or even complex ) arguments via the formula:

C_x = 4^x Gam(x+1/2) / [ sqrt(PI) Gam(x+2) ]

STACK	INPUTS	OUTPUTS
Level 1	x	C_x

Example:

PI CAT -> C_PI = 5.75308575086

-Likewise Cat (1+2.i) = 0.0452586883297 + 0.518182561743 i

i) Fresnel Integrals

-Fresnel integrals are defined by: S(x) = §₀^x sin(pi.t²/2).dt and C(x) = §₀^x cos(pi.t²/2).dt

Formulae:

S(x) = SUM_{n=0,1,2,.....} (-1)ⁿ (pi/2)²ⁿ⁺¹ x⁴ⁿ⁺³ / ((2n+1)!(4n+3))

C(x) = SUM_{n=0,1,2,.....} (-1)ⁿ (pi/2)²ⁿ x⁴ⁿ⁺¹ / ((2n)!(4n+1))

-A complex continued fraction is also employed if x > sqrt(Pi)

C(x) + i.S(x) = ((1+i)/2) erf z with z = (1-i).(x/2).(pi)^1/2

and ( 1 - erf z ) exp z² = (pi)^-1/2 ( 1/(z+0.5/(z+1/(z+1.5/(z+2/(z+ .... ))))) )

STACK	INPUTS	OUTPUTS
Level 2	/	C(x)
Level 1	x	S(x)

Examples:

C(1.5) = 0.445261176033
1.5 CSX -> S(1.5) = 0.697504960084

C(4) = 0.498426033044
4    CSX -> S(4) = 0.420515754248

     j) Debye Function

Definition: db(x;n) = §_x^+infinity tⁿ/(e^t-1).dt      where n is a positive integer and x > 0

Formula: db(x;n) = Sum_k>0 e^-k.x [ xⁿ/k + n.x^n-1/k² + ..... + n!/kⁿ⁺¹ ]

STACK	INPUTS	OUTPUTS
Level 2	x	/
Level 1	n	db(x;n)

Example:

0.7 ENTER
3 DBXN -> db( 0.7 ; 3 ) = 6.40683359558

k) Parabolic Cylinder Function

Formula:

D_n(x) = 2^n/2 Pi^1/2 exp(-x²/4) [ 1/Gam((1-n)/2) M( -n/2 , 1/2 , x²/2 ) - 2^1/2 ( x / Gam(-n/2) ) M [ (1-n)/2 , 3/2 , x²/2 ]

STACK	INPUTS	OUTPUTS
Level 2	n	/
Level 1	x	D_n(x)

Example:

   0.4   ENTER
   1.8   DNX    -> D_0.4(1.8) = 0.579579486208

     l) Error Function

Definition:

   erf x = (2/pi^1/2) §₀^x e^{-t^2} dt

Formula:

   erf x = (2x/Pi^1/2) exp(-x²) M( 1 , 3/2 , x² )

STACK	INPUT	OUTPUT
Level 1	x	erf(x)

Example:

2 ERF -> erf (2) = 0.995322265031

m) Fibonacci Functions

F_n(x) = [ 2^-n ( x + sqrt( x² +4 ) )ⁿ - Cos ( n.Pi ) 2ⁿ ( x + sqrt( x² +4 ) )^-n ] / sqrt ( x² +4 )

STACK	INPUTS	OUTPUTS
Level 2	n	/
Level 1	x	F_n(x)

Example:

1.2 ENTER
3.7 FIB -> F_1.2(3.7) = 1.27417856357

     n) Struve Functions

-The Struve function is defined by    H_n(x) = (x/2)ⁿ⁺¹ ₁F^~₂( 1 ; 3/2 , n + 3/2 ; - x²/4 )           where ₁F^~₂ is a regularized hypergeometric function.
and the modified Struve function:    L_n(x) = (x/2)ⁿ⁺¹ ₁F^~₂( 1 ; 3/2 , n + 3/2 ; x²/4 )

STACK	INPUTS	OUTPUTS
Level 2	n	/
Level 1	x	H_n(x) or L_n(x)

Example:

1.2 ENTER
3.4 HNX -> H_1.2 ( 3.4 ) = 1.11337265756

1.2 ENTER
3.4 LNX -> L_1.2 ( 3.4 ) = 4.64912946473

o) Lommel Functions

-"LOM1" & "LOM2" calculate the Lommel functions of the 1st & 2nd kind.
-They use the formulae:

s⁽¹⁾_m,n(x) = x^m+1 / [ (m+1)² - n² ] ₁F₂ ( 1 ; (m-n+3)/2 , (m+n+3)/2 ; -x²/4 )

    s⁽²⁾_m,n(x) = x^m+1 / [ (m+1)² - n² ] ₁F₂ ( 1 ; (m-n+3)/2 , (m+n+3)/2 ; x²/4 )
                      + 2^m+n-1 Gam(n) Gam((m+n+1)/2) x^-n / Gam((-m+n+1)/2) ₀F₁ ( ; 1-n ; -x²/4 )
                      + 2^m-n-1 Gam(-n) Gam((m-n+1)/2) xⁿ / Gam((-m-n+1)/2) ₀F₁ ( ; 1+n ; -x²/4 )

where _pF_q is the generalized hypergeometric function and Gam is Euler Gamma function.

STACK	INPUTS	OUTPUTS
Level 3	m	/
Level 2	n	/
Level 1	x	s⁽¹⁾_m,n(x) or s⁽²⁾_m,n(x)

Example:

   2   SQRT
   3   SQRT
PI LOM1 -> s⁽¹⁾_{sqrt(2),sqrt(3)}(PI) = 3.0030603835

   2   SQRT
   3   SQRT
PI LOM2 ->    s⁽²⁾_{sqrt(2),sqrt(3)}(PI) = 9.04879866075

     p) Polylogarithm

Definition:    Li_n(x) = SUM_k=1,2,... x^k/kⁿ

STACK	INPUTS	OUTPUTS
Level 2	n	/
Level 1	x	Li_n(x)

Example:

   3 ENTER
0.7 LIN -> Li₃(0.7) = 0.780063934258

     q) Incomplete Beta Function

-The incomplete Beta function is defined by   B_x(a,b) = §₀^x t^a-1 (1-t)^b-1 dt       ( a , b > 0 )
-We have also B_x(a,b) = (x^a/a) (1-x)^b F(1,a+b,a+1,x)   where F is the hypergeometric function

STACK	INPUTS	OUTPUTS
Level 3	x	/
Level 2	a	/
Level 1	b	B_x(a,b)

Example:

0.7 ENTER
    PI ENTER
1   E^X IBF ->    B_0.7(Pi,e) = 0.0296230460336

     r) Incomplete Gamma Function

-The incomplete gamma function is defind by    igam(a,x) = §₀^x e^-t t^a-1 dt

IGF employs:    igam(a,x) = (x^a/a).M(a,a+1,-x)         ( M = Kummer's Function )

STACK	INPUTS	OUTPUTS
Level 2	a	/
Level 1	x	igam(a,x)

Example:

     1.2 ENTER
   1.7 IGAM    ->    igam(1.2,1.7) = 0.697290896826

     s) Lerch Transcendent Function

-Lerch transcendent function is defined by    PHI( x , s , a ) = SUM_k=0,1,2,....   x^k / ( a + k )^s

STACK	INPUTS	OUTPUTS
Level 3	x	/
Level 2	s	/
Level 1	a	PHI( x , s , a )

Example:

0.7 ENTER
    PI   ENTER
   0.6 LERCH   -> F( 0.7 ; PI ; 0.6 ) = 5.1706011314

-Likewise,    PHI( 0.3 + 0.4 i , 3 + 4 i , 1 + 2 i ) = 7.65815909509 - 1.51511436459 i

Notes:

-Here, we assume that a # 0 , -1 , -2 , .....
-Execution time tends to infinity as | x | tends to 1.

     t) Coulomb Wave Functions

-The regular Coulomb wave function F_L(n,r) and the irregular Coulomb wave function G_L(n,r) are 2 independant solutions of the differential equation:

d²y/dr² + [ 1 - 2n/r - L(L+1)/r² ] y = 0

Formula:

F_L(n,r) = 2^L r^L+1 [ 1/Gam(2L+2) ] | Gam(L+1+i.n) | exp [ -(Pi) n/2 - i r ] ₁F₁ (L+1-i.n;2L+2;2i.r) where M_k,µ = Whittaker's function of the 1st kind

STACK	INPUTS	OUTPUTS
Level 3	L	/
Level 2	n	/
Level 1	r	F_L(n,r)

Example:

   1.4   ENTER
-2.1   ENTER
   1.6   RCWF -> F_1.4( -2.1 , 1.6 ) = 0.779739873373

Notes:

-In fact, with the example above, the HP48 returns   0.779739873373 + 3.247... E-11 i
-Without roundoff-error, imaginary part = 0

-The irregular Coulomb wave function can be computed by:

G_L(n,r) = i F_L(n,r) + | Gam(L+1+i.n) | [ 1/Gam(L+1+i.n) ] exp [ (PI) n/2 + i L(PI)/2 ] W_i.n,L+1/2 (2i.r)

where W_k,µ = Whittaker's function of the 2nd kind

STACK	INPUTS	OUTPUTS
Level 3	L	/
Level 2	n	/
Level 1	r	G_L(n,r)

Example:

   1.4   ENTER
-2.1   ENTER
   1.6   ICWF -> G_1.4( -2.1 , 1.6 ) = -0.20636826825

Notes:

-With the example above, the HP48 returns   0.77973987338 + 0.000000000039 i
-Without roundoff-error, imaginary part = 0

     u) Whittaker's Functions

-Whittaker's function of the 1st kind M_k,m(x) is defined by

M_k,m(x) = x^m+1/2 e^-x/2 M(m-k+1/2,2m+1,x)

where M(a,b,x) = Kummer's function.

STACK	INPUTS	OUTPUTS
Level 3	k	/
Level 2	m	/
Level 1	x	M_k,m(x)

Example:

2 SQRT
3 SQRT
PI WHIM -> M_{sqrt(2),sqrt(3)}(PI) = 5.61242619637

-Whittaker's function of the 2nd kind W_k,m(x) is defined by

W_k,m(x) = x^m+1/2 e^-x/2 U(m-k+1/2,2m+1,x)

where U(a,b,x) = hypergeometricU function.

STACK	INPUTS	OUTPUTS
Level 3	k	/
Level 2	m	/
Level 1	x	W_k,m(x)

Example:

2   SQRT
3   SQRT
PI WHIW -> W_{sqrt(2),sqrt(3)}(PI) = 2.17759342004

     v) Toronto Functions

-Toronto function is defined by

   T(m,n,r) = exp(-r²) [ Gam((m+1)/2) / n! ] r^2n-m+1 M( (m+1)/2 , n+1 , r² )

where M(a,b,x) = Kummer's function

STACK	INPUTS	OUTPUTS
Level 3	m	/
Level 2	n	/
Level 1	r	T(m,n,r)

Example:

   2   SQRT
   3   SQRT
PI   TOR -> T( sqrt(2),sqrt(3),PI ) = 0.963524225381

19°) Spheroidal Wave Functions

-We want to solve the differential equation ( 1 - x² ) d²S/dx² - 2.x dS/dx + [ L_mn - c² x² - m²/( 1 - x² ) S ] = 0

-First, LMN finds the eigenvalues L_mn by solving the transcendental equation U₁(L_mn) + U₂(L_mn) = 0 where U₁ & U₂ are 2 continued fractions:

                                                 -b_r^m                      -b_r-2^m
   U₁(L_mn) = g_r^m - L_mn +    ---------------       ---------------   ..............
                                           g_r-2^m - L_mn +         g_r-4^m - L_mn +

                                                                                                                           r = n - m
                                     -b_r+2^m                    -b_r+4^m
   U₂(L_mn) =           ---------------       ----------------   ..............
                               g_r+2^m - L_mn +         g_r+4^m - L_mn +

   with      b_r^m = [ r.(r-1).(2m+r).(2m+r-1).c⁴ ] / [ (2m+2r-1)².(2m+2r+1).(2m+2r-3) ]
   and      g_r^m = (m+r).(m+r+1) + (c²/2).[ 1 - (4m²-1)/((2m+2r-1).(2m+2r+3)) ]

-In fact, LMN only employs 12 iterations ( stored in 'ITER' ) to calculate the continued fractions.
-So, you can change the number of iterations - store another integer in 'ITER' - to check the accuracy

Example:      m = 0.2 , n = 0.6 , c² = 1.7

-Store these 3 numbers in the corresponding variables 'M' 'N' 'C2'

LMN   -> L_mn = 2.24686665135 ( we get the same result with 14 stored in 'ITER' )

-This value is stored automatically in 'L'

-Then SWF computes S_mn(x) with | x | < 1 by the following formulae:

S_mn(x) = ( 1 - x² )^m/2 f(x) where f(x) = Sum_k=0,1,.... a_k x^k

with a₀ = P_n^m(0) = 2^m sqrt(PI) / [ Gam((1-m-n)/2) Gam((2-m+n)/2 ] Flammer's scheme
and a₁ = P'_n^m(0) = ( m + n ) 2^m sqrt(PI) / [ Gam((2-m-n)/2) Gam((1-m+n)/2 ]

-In both cases, (k+1)(k+2) a_k+2 - [ k ( k + 2m + 1 ) - L_mn + m ( m + 1 ) ] a_k - c² a_k-2 = 0

STACK	INPUTS	OUTPUTS
Level 1	x	S_mn(x)

Example:

0.7   SWF ->     S_mn(0.7) = 0.682645661102

20°) Ellipsoidal Wave Functions

-The ellipsoidal wave function W(x) satisfies the differential equation:

      W''(x) = W(x) [ a + b k² Sn² (x) + q k⁴ Sn⁴ (x) ]     where Sn (x) is a Jacobian elliptic function.

-EWF computes a solution of the form W(x) = 1 + A₁ t + A₂ t² + ............. + A_m t^m + ................

     where t = Sn² (x | k² )

-This leads to the recurrence relation

     ( 2m+1 ) (2m+2 ) A_m+1 = q k⁴ A_m-2 + A_m-1 k² [ b - 2 ( m-1 ) ( 2m - 1 ) ] + A_m [ -a + 4 m² ( k² + 1 ) ]

STACK	INPUT	OUTPUT
Level 1	x	W(x)

Example: a = -1.2 b = 4.59 k² = 0.8 q = - sqrt(2)

-Calculate W(1)

-Sore the 4 constants a = -1.2    b = 4.59 k² = 0.8   q = - sqrt(2) into the corresponding variables 'a' 'b' 'K2' 'q'

-Then 1   EWF -> W(1) = 0.640627307889

Note:

SNX computes the Jacobian elliptic function: sn ( u | m ) ( cf §17-b )

Example: 0.7 ENTER 0.3 SNX -> sn ( 0.7 | 0.3 ) = 0.63230477631

21°) Euler Gamma Constant

'EULER' simply contains Euler's gamma constant = 0.577215664902

This is called by several programs

22°) Asymptotic Expansions

     a) Airy Functions

AIBI computes Ai(x) & Bi(x) for large positive x-values:

Ai(x) ~ { 1/2sqrt[ pi sqrt(x) ] } exp ( -2/3 x^3/2 ) ₂F₀( 1/6 , 5/6 ; -3/(4x^3/2) )
Bi(x) ~ { 1/sqrt[ pi sqrt(x) ] } exp ( 2/3 x^3/2 ) ₂F₀( 1/6 , 5/6 ; 3/(4x^3/2) )

STACK	INPUTS	OUTPUTS
Level 2	/	Ai(x)
Level 1	x > 0	Bi(x)

Example:

Ai(8) = 4.69220761616 E-8
8 AIBI -> Bi(8) = 1199586.00411

Note:

-If x is not large enough, the ( divergent ) series leads to an infinite loop.
-Press ON to stop the loops.

b) Coulomb Wave Functions

Formulae:

F_L(n,r) = g cos µ + f sin µ where µ = r - n Ln (2r) - L.(Pi)/2 + Arg Gam ( 1 + L + i n )
G_L(n,r) = f cos µ - g sin µ and f + i.g ~ u₀ + u₁ + u₂ + ............ + u_k + ......................

with u₀ = 1 and u_k / u_k-1 = ( i n - L + k - 1 ) ( i n + L + k ) / ( 2k i r ) ( k > 0 )

STACK	INPUTS	OUTPUTS
Level 3	L	/
Level 2	n	/
Level 1	r	F_L(n,r)

Example: L = 3 n = 4 r = 21

3 ENTER
4 ENTER F₃(4;21) = 1.13300294999
21 CWF -> G₃(4;21) = 0.121760594375

Note:

-If x is not large enough, the ( divergent ) series leads to an infinite loop.
-Press ON to stop the loops.

c) Parabolic Cylinder Functions

xⁿ exp(-x²/4) [ 1 - n(n-1) / ( 2 x² ) + n(n-1)(n-2)(n-3) / ( 2 ( 4 x⁴ ) ) - ....... ]

STACK	INPUTS	OUTPUTS
Level 2	n	/
Level 1	x	D_n(x)

Example:

   PI    ENTER
5.37 DNX    -> D_Pi(5.37) = 0.12839987247

Note:

-If x is not large enough, the ( divergent ) series leads to an infinite loop.
-Press ON to stop the loops.

     d) Exponential Integrals En

Formula: E_n(x) ~ (1/x) exp(-x) ₂F₀ ( 1 , n ; -1/x )

STACK	INPUTS	OUTPUTS
Level 2	n	/
Level 1	x	E_n(x)

Example:

   PI    ENTER
41 ENX    -> E_Pi(41) = 3.54610150354 E-20

Note:

-If x is not large enough, the ( divergent ) series leads to an infinite loop.
-Press ON to stop the loops.

     e) Struve Functions Hn

-We have:      H_n(x) ~ Y_n(x) + [ 2^1-n x^n-1 / sqrt(Pi) / Gam(n+1/2) ] ₃F₀( 1/2 , -n + 1/2 , 1 ; -4/x² )

STACK	INPUTS	OUTPUTS
Level 2	n	/
Level 1	x	H_n(x)

Example:

PI ENTER
24 ENX -> H_Pi(24) = 29.8495813024

Note:

-If x is not large enough, the ( divergent ) series leads to an infinite loop.
-Press ON to stop the loops.

f) Bessel Functions J & Y and K

Formulae:

J_n(x) = (2/(pi*x))^1/2 ( P.cos(x-(2n+1)pi/4) - Q.sin(x-(2n+1)pi/4) )
Y_n(x) = (2/(pi*x))^1/2 ( P.sin(x-(2n+1)pi/4) + Q.cos(x-(2n+1)pi/4) )

where P ~ 1 - (4n²-1)(4n²-9)/(2!(8x)²) + (4n²-1)(4n²-9)(4n²-25)(4n²-49)/(4!(8x)⁴) - ......
and Q ~ (4n²-1)/(8x) - (4n²-1)(4n²-9)(4n²-25)/(3!(8x)³) + ......

STACK	INPUTS	OUTPUTS
Level 2	n	J_n(x)
Level 1	x	Y_n(x)

Example:

PI ENTER J_Pi(16) = -0.0811607476181
16 JYNX -> Y_Pi(16) = -0.184298980133

Notes:

-In fact, the continued fraction used in JYNX in the BESSEL directory also works for large arguments

-With the modified Bessel function K_n(x) , we have the formula:

K_n(x) ~ (pi/(2x))^1/2 e^-x ( 1 + (4n²-1)/(8x) + (4n²-1)(4n²-9)/(2!(8x)²) + (4n²-1)(4n²-9)(4n²-25)/(3!(8x)³) + ...... )

STACK	INPUTS	OUTPUTS
Level 2	n	/
Level 1	x	K_n(x)

Example:

1.4 ENTER
19 KNX -> K_1.4(19) = 1.683198846503 E-9

Note:

-If x is not large enough, the ( divergent ) series leads to an infinite loop.
-Press ON to stop the loops.

23°) Fractional-Integro Differentiation

-If a function f is defined by a power series: f(x) = SUM_{k=0,1,2,.....} c_k x^k , its fractional integro-differentiation may be computed by

d^µ f / dx^µ = D^µ f(x) = SUM_{k=0,1,2,.....} c_k [ Gam(k+1) / Gam(k+1-µ) ] x^k-µ where µ is a real number ( integer or fractional )

-If the function may be expressed in terms of hypergeometric functions _pF_q , the following relation is very useful too:

D^µ _pF_q ( a₁ , .......... , a_p ; b₁ , .......... b_q ; x ) = x^-µ Gam(b₁).........Gam(b_q) _p+1F^~_q+1 ( 1 , a₁ , .......... , a_p ; 1-µ , b₁ , .......... b_q ; x )

   where   Gam = Euler's Gamma function and _pF^~_q is the regularized generalized hypergeometric function

     a) Exponential

    D^µ Exp x = x^-µ ₁F^~₁ ( 1 ; 1-µ ; x )

STACK	INPUTS	OUTPUTS
Level 2	µ	/
Level 1	x	D^µ Exp x

Example:

3.14 ENTER
1.28 DEXP -> D^3.14 Exp ( 1.28 ) = 3.50196367142

b) Logarithm

D^µ Ln x = x^-µ FC^(µ)_log (x)

where FC^(µ)_log (x) = (-1)^µ-1 (µ-1) ! if µ is a positive integer
and FC^(µ)_log (x) = [ Ln x - Psi(1-µ) - gamma ] / Gam(1-µ) otherwise

Psi = Digamma Function , gamma = Euler's constant = 0.5772156649... and Gam = Gamma Function.

STACK	INPUTS	OUTPUTS
Level 2	µ	/
Level 1	x	D^µ Ln x

Example:

     3.14   ENTER
     1.28 DLN -> D^3.14 Ln ( 1.28 ) = 1.13856985083

     c) FClog

FC^(µ)_log (x) = (-1)^µ-1 (µ-1) !                                               if µ is a positive integer
          and       FC^(µ)_log (x) = [ Ln x - Psi(1-µ) - gamma ] / Gam(1-µ)        otherwise

Psi = Digamma Function , gamma = Euler's constant = 0.5772156649... and Gam = Gamma Function.

STACK	INPUTS	OUTPUTS
Level 2	µ	/
Level 1	x	FC^(µ)_log (x)

Example:

3.14 ENTER
1.28 FCLN -> FC^(3.14)_log (1.28) = 2.47171836412

d) Cosine

D^µ Cos x = (2/x)^µ sqrt(PI) ₁F^~₂ ( 1 ; (1-µ)/2 , (2-µ)/2 ; -x²/4 )

STACK	INPUTS	OUTPUTS
Level 2	µ	/
Level 1	x	D^µ Cos x

Example:

     3.14   ENTER
     1.28 DCOS -> D^3.14 Cos ( 1.28 ) = 0.888787267089

     e) Sine

    D^µ Sin x = 2^µ-1 sqrt(PI) x^1-µ₁F^~₂ ( 1 ; (2-µ)/2 , (3-µ)/2 ; -x²/4 )

STACK	INPUTS	OUTPUTS
Level 2	µ	/
Level 1	x	D^µ Sin x

Example:

3.14 ENTER
1.28 DSIN -> D^3.14 Sin ( 1.28 ) = -1.91420925734 E-2

f) Hyperbolic Cosine

D^µ Cosh x = (2/x)^µ sqrt(PI) ₁F^~₂ ( 1 ; (1-µ)/2 , (2-µ)/2 ; x²/4 )

STACK	INPUTS	OUTPUTS
Level 2	µ	/
Level 1	x	D^µ Cosh x

Example:

     3.14   ENTER
     1.28 DCH -> D^3.14 Cosh ( 1.28 ) = 1.50295822021

     g) Hyperbolic Sine

     D^µ Sinh x = 2^µ-1 sqrt(PI) x^1-µ₁F^~₂ ( 1 ; (2-µ)/2 , (3-µ)/2 ; x²/4 )

STACK	INPUTS	OUTPUTS
Level 2	µ	/
Level 1	x	D^µ Sinh x

Example:

     3.14   ENTER
     1.28 DSH -> D^3.14 Sinh ( 1.28 ) = 1.9990054512

     h) Airy Functions

D^µ Ai(x) = 3^µ-4/3 x^-µ { 3^2/3 Gam(1/3) ₂F^~₃ [ 1/3 , 1 ; (1-µ)/3 , (2-µ)/3 , (3-µ)/3 ; x³/9 ] - x Gam(2/3) ₂F^~₃ [ 2/3 , 1 ; (4-µ)/3 , (2-µ)/3 , (3-µ)/3 ; x³/9 ] }

D^µ Bi(x) = 3^µ-5/6 x^-µ { 3^2/3 Gam(1/3) ₂F^~₃ [ 1/3 , 1 ; (1-µ)/3 , (2-µ)/3 , (3-µ)/3 ; x³/9 ] + x Gam(2/3) ₂F^~₃ [ 2/3 , 1 ; (4-µ)/3 , (2-µ)/3 , (3-µ)/3 ; x³/9 ] }

STACK	INPUTS	OUTPUTS
Level 2	µ	D^µ Ai(x)
Level 1	x	D^µ Bi(x)

Example:

     3.14   ENTER D^3.14 Ai ( 1.28 ) = -0.162004855985
     1.28 DAIBI ->    D^3.14 Bi ( 1.28 ) = 3.43259262532

     i) Error Function

D^µ Erf (x) = 2^µ x^1-µ ₂F^~₂ [ 1/2 , 1 ; (2-µ)/2 , (3-µ)/2 ; -x² ]

STACK	INPUTS	OUTPUTS
Level 2	µ	/
Level 1	x	D^µ Erf(x)

Example:

     3.14   ENTER
     1.28 DERF -> D^3.14 Erf ( 1.28 ) = 1.25055702267

     j) Hermite Functions

     D^µ H_n (x) = [ 2^n+µ (PI) x^-µ / Gam((1-n)/2) ] ₂F^~₂ [ 1 , -n/2 ; (1-µ)/2 , (2-µ)/2 ; x² ] - [ 2^n+µ (PI) x^1-µ / Gam((-n)/2) ] ₂F^~₂ [ 1 , (1-n)/2 ; 1-µ/2 , (3-µ)/2 ; x² ]

STACK	INPUTS	OUTPUTS
Level 3	µ	/
Level 2	n	/
Level 1	x	D^µ H_n (x)

Example:

     3.14   ENTER
     2.41   ENTER
     1.28   DHMT ->    D^3.14 H_2.41 ( 1.28 ) = 3.53770762114

k) Kummer's Function

D^µ F(a;b;x) = x^-µ Gam(b) ₂F^~₂ ( 1 , a ; 1-µ , b ; x )

STACK	INPUTS	OUTPUTS
Level 4	µ	/
Level 2	a	/
Level 2	b	/
Level 1	x	D^µ F(a;b;x)

Example:

     3.14   ENTER
     2 SQRT
3 SQRT
     1.28   DKUM -> D^3.14 F ( 2^1/2 ; 3^1/2 ; 1.28 ) = 2.07589150484

     l) Generalized Laguerre Functions

D^µ L^a_n (x) = [ Gam(n+a+1)/Gam(n+1) ] x^-µ₂F^~₂ ( 1 , -n ; a+1 , 1-µ ; x )

STACK	INPUTS	OUTPUTS
Level 4	µ	/
Level 2	a	/
Level 2	n	/
Level 1	x	D^µ L^a_n (x)

Example:

     3.14   ENTER
     1.76   ENTER
     2.41   ENTER
     1.28   DLANX -> D^3.14 L^1.76_2.41 ( 1.28 ) = -1.767203465

m) Bessel Functions

• Bessel Function - 1st kind D^µ J_n (x) = 2^µ-2n sqrt(PI) x^n-µ Gam(n+1) ₂F^~₃ [ (n+1)/2 , (n+2)/2 ; (n+1-µ)/2 , (n+2-µ)/2 , n+1 ; -x²/4 ]

• Bessel Function - 2nd kind

D^µ Y_n (x) = 2^µ-2n (PI)^1/2 x^-µ-n csc(n.PI) { -16ⁿ Gam(1-n) ₂F^~₃ [ (1-n)/2 , (2-n)/2 ; (1-µ-n)/2 , (2-µ-n)/2 , 1-n ; -x²/4 ]

+ x²ⁿ Cos(n.PI) Gam(n+1) ₂F^~₃ [ (n+1)/2 , (n+2)/2 ; (n+1-µ)/2 , (n+2-µ)/2 , n+1 ; -x²/4 ] } where n is not an integer.

• Modified Bessel Function - 1st kind D^µ I_n (x) = 2^µ-2n sqrt(PI) x^n-µ Gam(n+1) ₂F^~₃ [ (n+1)/2 , (n+2)/2 ; (n+1-µ)/2 , (n+2-µ)/2 , n+1 ; x²/4 ]

• Modified Bessel Function - 2nd kind

D^µ K_n (x) = 2^µ-2n-1 (PI)^3/2 x^-µ-n csc(n.PI) { 16ⁿ Gam(1-n) ₂F^~₃ [ (1-n)/2 , (2-n)/2 ; (1-µ-n)/2 , (2-µ-n)/2 , 1-n ; x²/4 ]

- x²ⁿ Gam(n+1) ₂F^~₃ [ (n+1)/2 , (n+2)/2 ; (n+1-µ)/2 , (n+2-µ)/2 , n+1 ; x²/4 ] } where n is not an integer.

• Spherical Bessel Function - 1st kind D^µ j_n (x) = 2^µ-2n-1 PI x^n-µ Gam(n+1) ₂F^~₃ [ (n+1)/2 , (n+2)/2 ; (n+1-µ)/2 , (n+2-µ)/2 , n+3/2 ; -x²/4 ]

STACK	INPUTS	OUTPUTS
Level 3	µ	/
Level 2	n	/
Level 1	x	D^µ Bessel_n (x)

Examples:

     3.14   ENTER
     2.41   ENTER
     1.28 DJNX ->    D^3.14 J_2.41 ( 1.28 ) = -0.150524582429

     3.14   ENTER
     2.41   ENTER
     1.28 DYNX ->    D^3.14 Y_2.41 ( 1.28 ) = 25.493085779

     3.14   ENTER
     2.41   ENTER
     1.28 DINX ->       D^3.14 I_2.41 ( 1.28 ) = 0.352247278862

     3.14   ENTER
     2.41   ENTER
     1.28 DKNX ->    D^3.14 K_2.41 ( 1.28 ) = -38.9846931335

     3.14   ENTER
     2.41   ENTER
     1.28 DSB1 ->    D^3.14 j_2.41 ( 1.28 ) = -0.064451621889

     n) Fresnel Integrals

• Fresnel Cosine Integral D^µ C(x) = 2^2µ-3/2 PI^3/2 x^1-µ₃F^~₄ [ 1/4 , 3/4 , 1 ; (2-µ)/4 , (3-µ)/4 , (4-µ)/4 , (5-µ)/4 ; -(PI)² x⁴/16 ]

• Fresnel Sine Integral D^µ S(x) = 2^2µ-11/2 PI^5/2 x^3-µ₃F^~₄ [ 3/4 , 1 , 5/4 ; (4-µ)/4 , (5-µ)/4 , (6-µ)/4 , (6-µ)/4 ; -(PI)² x⁴/16 ]

STACK	INPUTS	OUTPUTS
Level 2	µ	/
Level 1	x	D^µ C(x) or D^µ S(x)

Examples:

3.14 ENTER
1.28 DCX -> D^3.14 C ( 1.28 ) = 16.9561224908

3.14 ENTER
1.28 DSX -> D^3.14 S ( 1.28 ) = -11.203027733

o) Exponential Integrals & Related Functions

• Exponential Integral D^µ Ei x = [ FC^(µ)_log (x) + gamma / Gam(1-µ) ] x^-µ + x^1-µ ₂F^~₂ ( 1 , 1 ; 2 , 2-µ ; x )

• Sine Integral D^µ Si x = 2^µ-2 PI x^1-µ₂F^~₃ ( 1/2 , 1 ; 3/2 , (2-µ)/2 , (3-µ)/2 ; -x²/4 )

• Hyperbolic Sine Integral D^µ Shi x = 2^µ-2 PI x^1-µ₂F^~₃ ( 1/2 , 1 ; 3/2 , (2-µ)/2 , (3-µ)/2 ; x²/4 )

• Cosine Integral D^µ Ci x = [ FC^(µ)_log (x) + gamma / Gam(1-µ) ] x^-µ - 2^µ-3 sqrt(PI) x^2-µ ₂F^~₃ ( 1 , 1 ; 2 , (3-µ)/2 , (4-µ)/2 ; -x²/4 )

• Hyperbolic Cosine Integral D^µ Chi x = [ FC^(µ)_log (x) + gamma / Gam(1-µ) ] x^-µ + 2^µ-3 sqrt(PI) x^2-µ ₂F^~₃ ( 1 , 1 ; 2 , (3-µ)/2 , (4-µ)/2 ; x²/4 )

STACK	INPUTS	OUTPUTS
Level 2	µ	/
Level 1	x	D^µ ... (x)

Examples:

3.14 ENTER
1.28 DEI -> D^3.14 Ei ( 1.28 ) = 1.98213560559

3.14 ENTER
1.28 DSI -> D^3.14 Si ( 1.28 ) = -0.045395643837

3.14 ENTER
1.28 DSHI -> D^3.14 Shi ( 1.28 ) = 0.57649521068

3.14 ENTER
1.28 DCI -> D^3.14 Ci ( 1.28 ) = 1.36732389602

3.14 ENTER
1.28 DCHI -> D^3.14 Chi ( 1.28 ) = 1.4056403949

Notes:

-The results are accurate only if x is enough small.
-S1 S2 ....... S7 are subroutines called by the programs above.

24°) Non-Linear Systems

     a) Secant Method

SOLVE employs the generalized secant method to solve a system of n equations in n unkowns.
They must be written under the form:

f₁( x₁,... ,x_n ) = 0
   .....................

f_n( x₁, ... ,x_n ) = 0

SOLVE takes a program in level 3 that takes the vector V = [ x₁ , ...... , x_n ] and returns   [ f₁(x_{1, ....,}x_n) , ...... , f_n(x_{1,
....,}x_n) ] in level 1

and 2 initial guesses [ x₁⁽¹⁾ , ...... , x_n⁽¹⁾ ]
and [ x₁⁽²⁾ , ...... , x_n⁽²⁾ ] in levels 1 & 2

-The successive approximations are displayed and - if the method is convergent - gives | F(V) | in level 2 & a solution in level 1 ( level 3 is unchanged )

STACK	INPUTS	OUTPUTS
Level 3	<< F >>	<< F >>
Level 2	[ x₁⁽¹⁾ , ...... , x_n⁽¹⁾ ]	\| F(V) \|
Level 1	[ x₁⁽²⁾ , ...... , x_n⁽²⁾ ]	V_Solution

Example: Find a solution near ( 4 , 1 , 3 , 6 ) of the system:

       x + y + z + t - 16 = 0
          x y z - 3 t = 0
       4 x² - y z t - 40    = 0
          x y z t - 140 = 0

<< OBJ-> DROP -> X Y Z T
<< X Y Z T + + + 16 -
X Y Z * * T 3 * -
X SQ 4 * Y Z T * * - 40 -
X Y Z T * * * 140 -
4 ->ARRY
>>
>> ENTER
[ 4 1 3 6 ] ENTER
[ 4.1 1.1 3.1 6.1 ] SOLVE the successive approximations are displayed and eventually, we get:

Level 2 0
Level 1 [ 4.26654047494 1.35363223606 3.54852677842 6.83130051063 ]

-So f(V) ~ 0 and we have one solution:

x = 4.26654047494
y = 1.35363223606
z = 3.54852677842
t = 6.83130051063

Note:

-The iterations stop when 2 successive x-values are identical, so
-All the coordinates of the 2 initial guesses must be different: x_j⁽¹⁾ # x_j⁽²⁾ for all j

b) Newton's Method

NEWT also solves systems:

f₁( x₁,... ,x_n ) = 0
.....................

f_n( x₁, ... ,x_n ) = 0

NEWT takes the list of n functions in level 3
the list of n variables in level 2
and the initial guess in level 1

-First, it computes all the partial derivatives and then employs Newton's method to find a solution

STACK	INPUTS	OUTPUTS
Level 3	{ 'f₁' ........... 'f₁' }	{ 'f_j' + partial deriv }
Level 2	{ 'x₁' ......... 'x₁' }	{ 'x₁' ......... 'x₁' }
Level 1	[ x₁⁽²⁾ , ...... , x_n⁽²⁾ ]	V_Solution

Example: Find a solution near ( 4 , 1 , 3 , 6 ) of the system:

       x + y + z + t - 16 = 0
          x y z - 3 t = 0
       4 x² - y z t - 40    = 0
          x y z t - 140 = 0

{ 'X+Y+Z+T-16' 'X*Y*Z-3*T' '4*X^2-Y*Z*T-40' 'X*Y*Z*T-140' } ENTER
{ X Y Z T } ENTER
[ 4 1 3 6 ] NEWT

-After displaying the successive approximations, we get in level 1 [ 4.26654047494 1.35363223606 3.54852677842 6.83130051063 ]

-So we have one solution:

x = 4.26654047494
y = 1.35363223610
z = 3.54852677833
t = 6.83130051063

Notes:

-The list in level 3 now contains also all the partial derivatives.
-So, if you want to find another solution of the same system, don't modify level 3:
-Simply replace the vector in level 1 by another guess.

-The variable 'eps' contains a very small number ( 10^-10 ) to determine the end of the loops.

25°) Continued Fractions

CFR computes a function f(x) defined by a continued fraction:

                            a₁       a₂       a₃                               a_n
      f(x) = b₀ + -----   -----   -----   .....................   ----- .......................     where a_n & b_n are functions of x and n
                          b₁ +    b₂ +    b₃ +                            b_n +

-Themodified Lenz's algorithm is employed:

    we set   f₀ = b₀ if b₀ # 0 or f₀ = tiny if b₀ = 0    ( tiny = E-41 in CFR )
               C₀ = f₀ and D₀ = 0

then, for n = 1 , 2 , ....... C_n = b_n + a_n/C_n-1 ( C_n is replaced by tiny if C_n = 0 )
D_n = 1/ [ b_n + a_n D_n-1 ] ( D_n = 1/tiny if the denominator = 0 )

and f_n = C_n D_n f_n-1 until | C_n D_n - 1 | <= a "small" number ( E-11 in CFR )

>>>> f(x) = the last f_n

STACK	INPUTS	OUTPUTS
Level3	b₀	/
Level 2	<< F >>	/
Level 1	x	f(x)

Where b₀ is the 1st term of the continued fraction

and << F >> is a program that takes x in level 2 & n in level 1 and returns a_n in level 2 & b_n in leve 1 without changing the other stack levels

Example:    f is defined by:    b₀ = 1 ;   a_n = 2.x + n , b_n = x² + n² ( n > 0 ) ;    evaluate f(Pi)

   1 ENTER
<< OVER 2 * OVER + ROT SQ ROT SQ + >>    ENTER
PI CFR -> f(Pi) = 1.63626566028

References:

[1] Abramowitz and Stegun - "Handbook of Mathematical Functions" - Dover Publications - ISBN 0-486-61272-4
[2] http://functions.wolfram.com/
[3] Press , Vetterling , Teukolsky , Flannery - "Numerical Recipes in C++" - Cambridge University Press - ISBN 0-521-75033-4
[4] P. Borwein - "An Efficient Algorithm for the Riemann Zeta Function"
[5] B.C. Carlson, "Numerical Computation of Real or Complex Elliptic Integrals"