-JMB MATH  Manual

OVERVIEW
 

-The first 15 functions are actually M-Code routines.
-They work like built-in functions, saving X in L-register.
-However, most of them do not check for overflow.

-The French names are used for the hyperbolic functions:   SH  CH  TH  ASH  ACH  ATH
 
 

XROM	Function	Desciption
04,00  04,01  04,02  04,03  04,04  04,05  04,06  04,07  04,08  04,09  04,10  04,11  04,12  04,13  04,14  04,15  04,16  04,17  04,18  04,19  04,20  04,21  04,22  04,23  04,24  04,25  04,26  04,27  04,28  04,29  04,30  04,31  04,32  04,33  04,34  04,35  04,36  04,37  04,38  04,39  04,40  04,41  04,42  04,43  04,44  04,45  04,46  04,47	-JMB MATH  SH  CH  TH  ASH  ACH  ATH  HGF+  1/GM  Y^X  X#Y??  RAND  -ELLIPTIC  RF  RFZ  RJ  RJZ  "RG"  "LEI1"  "LEI2"  "LEI3"  "JEF"  "SAE"  -ORTHOPOL  "LEG"  "HMT"  "CHB"  "CHB2"  "LANX"  "PMN"  "PABNX"  "USP"  -MISCELLA  "SIMPLEX"  "NCSI"  "3DLS"  "LAGR"  "LAGRI"  "BVI"  "SNK"  "CGW3J"  "W6J"  "W9J"  "DTI"  -THE EARTH  "TGD"  "GRV"  GEU	Section Header  Hyperbolic Sine  Hyperbolic Cosine  Hyperbolic Tangent  Inverse Hyperbolic Sine  Inverse Hyperbolic Cosine  Inverse Hyperbolic Tangent  Generalized Hypergeometric Functions  Reciprocal Gamma Function  Numerical Power  Similar Comparison  Random Number Generator (with Timer)  Section Header  Carlson Elliptic Integral 1st Kind - Real Arguments  Carlson Elliptic Integral 1st Kind - Complex Arguments  Carlson Elliptic Integral 3rd Kind - Real Arguments  Carlson Elliptic Integrals 3rd Kind - Complex Arguments  Carlson Symmetric Eliiptic Integral - 2nd Kind  Legendre Elliptic Integral 1st Kind  Legendre Elliptic Integral 2nd Kind  Legendre Elliptic Integral 3rd Kind  Jacobian Elliptic Functions   Surface Area of an Ellipsoid   Section Header  Legendre Polynomials  Hermite Polynomials  Chebyshev Polynomials - 1st Kind  Chebyshev Polynomials - 2nd Kind  Generalized Laguerre Polynomials  Associated Legendre Functions 1st Kind  Jacobi Polynomials  Ultraspherical Polynomials  Section Header  Simplex Method  Natural Cubic Spline Integration  Tri-Diagonal Systems  Lagrange Interpolation Formula  Integration of the Lagrange Polynomial  Bivariate Interpolation  Stirling Numbers - 1st & 2nd Kind  Clebsh-Gordan Coeffs. & Wigner 3j-symbol  Wigner 6j-symbol  Wigner 9j-symbol  Double & Triple Integrals  Section Header  Terrestrial Geodesic Distance  Ellipsoidal Gravity  Euler's Gamma Constant
05,00  05,01  05,02  05,03  05,04  05,05  05,06  05,07  05,08  05,09  05,10  05,11  05,12  05,13  05,14  05,15  05,16  05,17  05,18  05,19  05,20  05,21  05,22  05,23  05,24  05,25  05,26  05,27  05,28  05,29  05,30  05,31  05,32  05,33  05,34  05,35  05,36  05,37  05,38  05,39  05,40  05,41  05,42  05,43  05,44  05,45	-SPEC FCNS  "AIRY"  "ERF"  "CSX"  "SI"  "SHI"  "CI"  "CHI"  "CSIX>3"  "EI"  "ENX"  "ZETA"  "INX"  "JNX"  "JNX1"  "YNX"  "KNX"  "DNX"  "KLV"  "KLV2"  "HNX"  "LNX"  "BETA"  "BMN"  "IBF"  "IGF"  "WEBAN"  "DBY"  "ALF"  "ALF2"  "FDF"  "PX2N"  "STD"  "RCWF"  "TMNR"  "LMNC2"  "SWF"  "HGF"  "RHGF"  -ASYMPTOT!  "AEAIRY"  "AEJYNX"  "AEKNX"  "AEDNX"  "AEKLV"  "AEHNX"	Section Header  Airy Functions  Ai(x) & Bi(x)  Error Function   Fresnel Integrals   C(x) & S(x)   Sine Integral   Si(x)   Hyperbolic Sine Integral   Cosine Integral   Ci(x)   Hyperbolic Cosine Integral   Continued fraction for  Ci(x) & Si(x)  if x > 3  Exponential Integral   Ei(x)  Exponential Integral   En(x)  Riemann Zeta Function  Modified Bessel Function 1st Kind   Bessel Function 1st Kind   Bessel Function 1st Kind , Integer order , Miller Alg.  Bessel Function 2nd Kind   Modified Bessel Function 2nd Kind   Parabolic Cylinder  Dn(x)   Kelvin Functions 1st Kind    bern(x) & bein(x)  Kelvin Functions 2nd Kind   kern(x) & kein(x)  Struve Function H   Struve Function L  Beta function   Beta function - Integer Arguments  Incomplete beta function,   Bx(a,b)   Incomplete gamma function, gam(a,x)  Weber & Anger functions   Debye Functions   Associated Legendre function 1st Kind - Pmn(x)  Associated Legendre function 2nd Kind - Qmn(x)  F-Distribution   Q( x | n , m ) = UTPF ( n , m , x )   Chi2 Probability Function    P( X2 | m )   Student's t-Distribution   A( t | n ) & UTPT( t | n )   Regular Coulomb Wave Function   Toronto functions   T(m,n,r)   Spheroidal Eigenvalues  Angular Spheroidal Wave Function of 1st Kind  Hypergeometric function   Regularized hypergeometric function   Section Header  Airy Functions - Asymptotic Expansion   Bessel Functions 1st & 2nd Kind - Asymptotic Expan.  Modified Bessel Functions 2nd Kind - Asympt. Expan.  Parabolic Cylinder - Asymptotic Expansion  Kelvin Functions 1st & 2nd Kind - Asympt. Expansion   Struve Function H - Asymptotic Expansion

 

HYPERBOLIC SINE
 
 

STACK	INPUT	OUTPUT
X	x	Sinh(x)
L	/	x

  1.2   XEQ "SH"  >>>   Sinh(1.2) = 1.509461356

Note:

  241  XEQ "SH"  returns the - apparently - wrong result  2.311746150 E-96

  but if you press SQRT it yields 1.520442748 E52 , which is correct !

-In fact,  Sinh 241 = 2.311746150 E104  is properly calculated and handled by the HP-41 despite the misleading display.

-Likewise, Sinh 460 is correctly computed.
-However, Sinh 461 is not !
-So, avoiding the check for overflow may be useful, provided all the intermediate results remain < E200.
 

HYPERBOLIC COSINE
 
 

STACK	INPUT	OUTPUT
X	x	Cosh(x)
L	/	x

  1.2   XEQ "CH"  >>>   Cosh(1.2) = 1.810655568
 

HYPERBOLIC TANGENT
 
 

STACK	INPUT	OUTPUT
X	x	Tanh(x)
L	/	x

    1.2   XEQ "TH"  >>>   Tanh(1.2) = 0.833654607
 

INVERSE HYPERBOLIC SINE
 
 

STACK	INPUT	OUTPUT
X	x	Asinh(x)
L	/	x

    1.2   XEQ "ASH"  >>>   Asinh(1.2) = 1.015973134
 

INVERSE HYPERBOLIC COSINE
 
 

STACK	INPUT	OUTPUT
X	x	Acosh(x)
L	/	x

  1.2   XEQ "ACH"  >>>   Acosh(1.2) = 0.622362504
 

INVERSE HYPERBOLIC TANGENT
 
 

STACK	INPUT	OUTPUT
X	x	Atanh(x)
L	/	x

    0.7   XEQ "ATH"  >>>   Atanh(0.7) = 0.867300528
 

GENERALIZED HYPERGEOMETRIC FUNCTIONS
 

 HGF+  computes    _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   ( Pochhammer symbol )
                                   ~
                and           _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!
 
 

STACK	INPUTS	OUTPUTS
T	T	T
Z	m	last k-value
Y	+/- p	1st neglected term
X	x	f(x)
L	/	x

  where  f(x) =   _mF_p( a₁,a₂,....,a_m ; b₁,b₂,....,b_p ; x )    if   Y = + p      hypergeometric function
                  _~
   and     f(x) =   _mF_p( a₁,a₂,....,a_m ; b₁,b₂,....,b_p ; x )    if   Y = - p     regularized hypergeometric function

 >>>   a₁,a₂,....,a_m   are to be stored in registers  R01 ............ Rmm
 and    b₁,b₂,....,b_p   are to be stored in registers  Rm+1 ........ Rm+p

Examples:

    •    1  STO 01   4  STO 02   7  STO 03   2  STO 04   3  STO 05   6  STO 06   5  STO 07

    3   ENTER^
    4   ENTER^
   PI  XEQ "HGF+"  >>>>   ₃F₄( 1 , 4 , 7 ; 2 , 3 , 6 , 5 ; Pi )  =  1.631019643

    3   ENTER^
  -4   ENTER^                        ~
   PI  XEQ "HGF+"   >>>>   ₃F₄( 1 , 4 , 7 ; 2 , 3 , 6 , 5 ; Pi )  =  0.0002831631326

   •    1  STO 01   4  STO 02   -2  STO 03   -5  STO 04

   2    ENTER^
  -2   ENTER^                        ~
  0.1  XEQ "HGF+"   >>>>   ₂F₂( 1 , 4 ; -2 , -5 ; 0.1 ) = 0.01289656888

Notes:

-If m = p = 0 , HGF+  returns exp(x)
-The program doesn't check if the series are convergent or not.
-Even when they are convergent, execution time may be prohibitive: press any key to stop HGF+

-Stack register T is saved and x is saved in L-register.
-R00 is unused.
-The alpha "register" is cleared.

-This M-Code routine is extensively used to compute special functions...

Warning:

-This routine checks that Rm+n exists  but  HGF+ does not check for alpha data...
 

RECIPROCAL OF THE GAMMA FUNCTION
 

-The following asymptotic expansion is used:

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

  with the relation:  Gam(x+n) = (x+n-1) (x+n-2) ...... (x+1) x Gam(x)   until  x+n > 5
 
 

STACK	INPUT	OUTPUT
X	x	1/Gam(x)
L	/	x

  PI   XEQ "1/GM"  >>>  1/Gam(Pi) = 0.437055717
 

POWER FUNCTION
 

-The M-Code routine  Y^X  gives the same results as the built-in Y^X except that it returns  0^0 = 1
 

A NUMERICAL TEST:  X#Y??
 

-In a program, the following instruction is skipped if X is "almost" equal to Y:
-Useful to avoid infinite loops ( in Newton's method ... and so on ... )
-However, if X or Y = 0, the loop will stop iff  X = Y = 0
 

RANDOM NUMBER GENERATOR
 
 

STACK	INPUT	OUTPUT
X	x	0 < r < 1
L	/	x

-The M-Code routine "RAND" generates pseudo-random numbers in X.
-It uses the internal clock and saves X in L-register.
 

-ELLIPTIC
 

CARLSON ELLIPTIC INTEGRAL 1ST KIND, REAL ARGUMENTS
 

    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
 
 

STACK	INPUTS	OUTPUTS
Z	z	/
Y	y	/
X	x	RF(x,y,z)
L	/	x

    2   ENTER^
    3   ENTER^
    4   XEQ "RF"  >>>   R_F(4,3,2) = 0.584082842
 

CARLSON ELLIPTIC INTEGRAL 1ST KIND, CONJUGATE COMPLEX ARGUMENTS
 

  "RFZ"  calculates  R_F(x,y+i.z,y-i.z)
 
 

STACK	INPUTS	OUTPUTS
Z	z	/
Y	y	/
X	x	RF(x,y+i.z,y-i.z)
L	/	x

  2   ENTER^
  3   ENTER^
  4   XEQ "RFZ"  >>>   R_F(4,3+2.i,3-2.i) = 0.529597761
 

CARLSON ELLIPTIC INTEGRAL 3RD KIND, REAL ARGUMENTS
 

    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
 
 

STACK	INPUT	OUTPUT
T	p > 0	/
Z	z	/
Y	y	/
X	x	RJ(x,y,z,p)

  1   ENTER^
  2   ENTER^
  3   ENTER^
  4   XEQ "RJ"  >>>   R_J(4,3,2,1) = 0.360378094

-Alpha register is cleared.
 

CARLSON ELLIPTIC INTEGRAL 3RD KIND, CONJUGATE COMPLEX ARGUMENTS
 

    "RJZ" computes  R_J(x,y+i.z,y-i.z,p)  with p > 0
 
 

STACK	INPUTS	OUTPUT
T	p > 0	/
Z	z	/
Y	y	/
X	x	RJ(x,y+i.z,y-i.z,p)

  1   ENTER^
  2   ENTER^
  3   ENTER^
  4   XEQ "RJZ"  >>>   R_J(4,3+2.i,3-2.i,1) = 0.287311651

-Alpha register is cleared.
 

  ·  All the following functions are focal programs - except Euler's constant
 
 

CARLSON SYMMETRIC ELLIPTIC INTEGRAL 2ND KIND, REAL ARGUMENTS
 

     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
 

REGISTERS:    R00 thru R03
FLAGS:            /
 
 

STACK	INPUTS	OUTPUT
Z	z	/
Y	y	/
X	x	RG(x,y,z)

  2   ENTER^
  3   ENTER^
  4   XEQ "RG"  >>>   R_G(4,3,2) = 1.725503029
 

LEGENDRE ELLIPTIC INTEGRAL 1ST KIND
 

    F ( x | m ) =  §₀^x ( 1 - m sin² t ) ^-1/2 dt  = sin x . R_F ( cos² x ; 1 - m sin² x ; 1 )
 

REGISTERS:    R00 thru R03
FLAGS:            F01-F02
 
 

STACK	INPUTS	OUTPUT
Y	m	/
X	x	F ( x | m )

  >>>   x is expressed in the current angular mode

   in DEG mode

   0.7   ENTER^
   84°  XEQ "LEI1"  >>>   F( 84° | 0.7 ) = 1.884976270

Note:

  x must be between  -90° ( exclusive ) and +90° ( inclusive )
 

LEGENDRE ELLIPTIC INTEGRAL 2ND KIND
 

     E ( x | m ) =  §₀^x ( 1 - m sin² t ) ^1/2 dt  = sin x . R_F ( cos² x ; 1 - m.sin² x ; 1 ) - (m/3) sin³ x . R_J ( cos² x ; 1 - m.sin² x ; 1 ; 1 )
 

REGISTERS:    R00 thru R04
FLAGS:            F01-F02
 
 

STACK	INPUTS	OUTPUT
Y	m	/
X	x	E ( x | m )

  >>>   x is expressed in the current angular mode

   in DEG mode

   0.7   ENTER^
   84°  XEQ "LEI2"  >>>   1.184070047

Note:

  x must be between  -90° ( exclusive ) and +90° ( inclusive )
 

LEGENDRE ELLIPTIC INTEGRAL 3RD KIND
 

     ¶ ( n ; x | m ) =  §₀^x  ( 1 + n sin² t ) ^-1 ( 1 - m sin² t ) ^-1/2 dt
                         = sin x . R_F ( cos² x ; 1 - m.sin² x ; 1 ) - (n/3) sin³ x . R_J ( cos² x ; 1 - m.sin² x ; 1 ; 1 + n.sin² x )
 

REGISTERS:    R00 thru R04
FLAGS:            F01-F02
 
 

STACK	INPUTS	OUTPUT
Z	n	/
Y	m	/
X	x	¶ ( n ; x | m )

  >>>   x is expressed in the current angular mode

   in DEG mode

   0.9   ENTER^
   0.7   ENTER^
   84°  XEQ "LEI3"  >>>   1.336853615

Note:

  x must be between  -90° ( exclusive ) and +90° ( inclusive )
-The sign convention for n is opposite that of Abramowitz & Stegun
 

JACOBIAN ELLIPTIC FUNCTIONS
 

 Let  m  be a number that verifies:  0 < = m < = 1      ( m is called the parameter )

   if    u =  §₀^v ( 1 - m sin² t ) ^-1/2 dt     ( 0 < = v < = 90° )    the angle  v  = am u    is called the amplitude,

 >>>   The 3 Jacobian elliptic functions  sn ; cn ; dn  are defined by:

      sn ( u | m ) = sin v  ; cn ( u | m) = cos v ;  dn ( u | m ) = ( 1 - m sin² v )^1/2

-The 9 other Jacobian elliptic functions can be obtained with the relations:

    cd = cn / dn  ;  dc = dn / cn  ;  ns = 1 / sn
    sd = sn / dn  ;   nc = 1 / cn   ;  ds = dn / sn
    nd = 1 / dn   ;   sc = sn / cn  ;  cs = cn / sn
 

REGISTERS:    R00 thru R07
FLAGS:            F06-F07
 
 

STACK	INPUTS	OUTPUTS
Z	/	dn ( u | m )
Y	m	cn ( u | m )
X	u	sn ( u | m )

Examples:

       0.3 ENTER^
       0.7 XEQ "JEF"    gives        sn ( 0.7 | 0.3 )  = 0.632304777
                                  RDN        cn ( 0.7 | 0.3 ) = 0.774719736
                                  RDN        dn ( 0.7 | 0.3 ) = 0.938113640

    Likewise,                 sn ( 0.7 | 1 ) = 0.604367777          sn ( 0.7 | 2 ) = 0.564297008        sn ( 0.7 | -3 ) = 0.759113422
                                   cn ( 0.7 | 1 ) = 0.796705460          cn ( 0.7 | 2 ) = 0.825571855        cn ( 0.7 | -3 ) = 0.650958381
                                   dn ( 0.7 | 1 ) = 0.796705460          dn ( 0.7 | 2 ) = 0.602609139        dn ( 0.7 | -3 ) =1.651895747
 

SURFACE AREA OF AN ELLIPSOID
 

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

REGISTERS:    R00 thru R03
FLAGS:            /
 
 

STACK	INPUTS	OUTPUTS
Z	c	/
Y	b	/
X	a	Area

Example:

    2   ENTER^
    4   ENTER^
    7   XEQ "SAE"  >>>  224.5830064
 

-ORTHOPOL
 

-In the following programs that computes orthogonal polynomials, n must be a positive integer smaller than 1000
 

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
 

REGISTERS:    R00 thru R02
FLAGS:            /
 
 

STACK	INPUTS	OUTPUTS
Y	n	Pn-1(x)
X	x	Pn(x)

       ( 0 < n < 1000 ,  n integer )

Example:

    7    ENTER^
  4.9   XEQ "LEG"  gives     P₇(4.9) =  1698444.019
                               RDN    P₆(4.9) = 188641.3852  is in Y-register.
 

HERMITE POLYNOMIALS
 

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

REGISTERS:    R00-R01
FLAGS:            /
 
 

STACK	INPUTS	OUTPUTS
Y	n	Hn-1(x)
X	x	Hn(x)

       ( 0 < n < 1000 ,  n integer )

Example:

   7       ENTER^
  3.14  XEQ "HMT"  >>>   H₇ (3.14) =   73726.24332
                                RDN   H₆ (3.14) =   21659.28040
 

CHEBYSHEV POLYNOMIALS - 1ST KIND
 

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

REGISTERS:    R00-R01
FLAGS:            F02
 
 

STACK	INPUTS	OUTPUTS
Y	n	Tn-1(x)
X	x	Tn(x)

       ( 0 < n < 1000 ,  n integer )

Example:

     7       ENTER^
  0.314  XEQ "CHB"  >>>>   T₇(0.314) =  -0.786900700
                                  RDN    T₆ (0.314) =    0.338782777
 

CHEBYSHEV POLYNOMIALS - 2ND KIND
 

       Formulae:    U_n(x) = 2x.U_n-1(x) - U_n-2(x)  ;  U₀(x) = 1  ;  U₁(x) = 2x
 

REGISTERS:    R00-R01
FLAGS:            F02
 
 

STACK	INPUTS	OUTPUTS
Y	n	Un-1(x)
X	x	Un(x)

       ( 0 < n < 1000 ,  n integer )

Example:

     7       ENTER^
  0.314   XEQ "CHB2"     >>>>    U₇ (0.314) =  -0.582815681
                                         RDN    U₆ (0.314) =   0.649952293
 

GENERALIZED LAGUERRE POLYNOMIALS
 

        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)
 

REGISTERS:    R00 thru R03
FLAGS:            /
 
 

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"  >>>>   L₇^(1.4)(Pi) = 1.688893513
                                 RDN    L₆^(1.4)(Pi) = 2.271353727
 

ASSOCIATED LEGENDRE FUNCTIONS 1st KIND - INTEGER ORDER & DEGREE
 

       Formulae:      (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
 

 REGISTERS:    R00 thru R03
 FLAGS:            /
 
 

STACK	INPUTS	OUTPUTS
Z	m	/
Y	n	Pnm-1(x)
X	x	Pnm(x)

    ( 0 £ m £ n )

Examples:

   4  ENTER^
   7  ENTER^
   3  XEQ "PMN"  >>>>     P₇⁴(3) =  37920960
                             X<>Y    P₆⁴(3) =  2963520

     3   ENTER^
  100  ENTER^
   0.7  XEQ "PMN"  >>>>   P₁₀₀³(0.7) = -58239.28386
                                X<>Y    P₉₉³(0.7) =  13003.91702

-Actually,  P_n^m(x)  is not always a polynomial !
 

JACOBI POLYNOMIALS
 

       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)
 

REGISTERS:    R00 thru R09
FLAGS:            /
 
 

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 "PABNX"  >>>>   P₇^(1.4;1.7)(1/Pi) = -0.322234421
                                           X<>Y  P₆^(1.4;1.7)(1/Pi) =   0.538220923
 

ULTRASPHERICAL POLYNOMIALS
 

           Formulae:      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)
 

REGISTERS:    R00 thru R03
FLAGS:            F02 ( if a = 0 )
 
 

STACK	INPUTS	OUTPUTS
Z	a	/
Y	n	Cn-1(a)(x)
X	x	Cn(a)(x)

       ( 0 < n < 1000 ,  n integer )

 C_n-1^(a)(x)  is in Y-register only if  a # 0

Examples:

  1.5  ENTER^
   7    PI   1/X   XEQ "USP"  >>>>   C₇^(1.5)(1/Pi)  =  -0.989046386
                                            X<>Y   C₆^(1.5)(1/Pi)  =    1.768780932

   0   ENTER^
   7   PI   1/X   XEQ "USP"  >>>>   C₇⁽⁰⁾(1/Pi)  =  -0.219109374
 

-MISCELLA
 

LINEAR PROGRAMMING: SIMPLEX METHOD
 

   SIMPLEX solves linear programming problems:

The purpose is to find  m  non-negative real numbers:   x₁ , ..... , x_m
satisfying:
                          b_i ³ a_i;1x₁ + ....... + a_i;m x_m            i = 1 , .... , n                                ( n inequations )                   all the b_i  , b_i'  and  b_i"
                          b_i' = a_i';1x₁ + ....... +a_i';m x_m            i' = 1 , ..... , n'                             ( n' equations )                    must be non-negative*
                          b_i"£ a_i";1x₁+ ....... +a_i";m x_m         i" = 1 , ..... , n"                             (n" inequations )

and such that      F  =  c₁x₁ + ....... +  c_mx_m       is maximized.    ( to find the minimum of F, seek the maximum of  -F )
 

 *  if some b_i are negative, multiply the corresponding lines by -1.

-Always write ( or re-write ) the system under the form above.
-Store the coefficients in column order from register R01 with 0 instead of F.
 

Warning:
 

-All the c_j  must be significantly smaller than 100.
  As an example if  F = 2400 x + 1200 y  it would be better to find the maximum of  2.4 x + 1.2 y and to multiply the result by 1000.

-This recommendation does not apply if there are only inequations of the first type ( ³ )  i-e if  n' = n" = 0.
-This program finds only 1 solution ( even if there are several solutions ).
-Status register "a"  is used.
-Therefore, "SIMPLEX" must not be called as more than a first level subroutine.
 

 REGISTERS:    SIZE = 1+(p+1)(m+n"+p+1 )     with p = n+n'+n"
 FLAGS:            F01-F02

   >>>  If F01 and/or F02 is set at the end, the problem has no solution.
 
 

STACK	INPUTS	OUTPUTS
T	n  = the number of inequations of the 1st kind ( ³ )	/
Z	n' = the number of equations ( = )	/
Y	n" = the number of inequations of the 2nd kind ( £ )	/
X	m = the number of unknowns	max (F)

Example:

Find  4 non-negative numbers  x, y, z and t satisfying:

             56 ³  x + 2y + 3z +4t
             41 ³ 5x +  y - 3z - t
             61 ³ 4x + 7y +2z - 3t
             25 = 2x + 3y - z + 2t                          ( always write the "³" inequations first, then the "=" equations , and finally the "£" inequations )
               7 £   x + y + z + t
             17 £ 2x - y + z + 3t

such that F = 2x + 4y +3z + 5t is maximized.
 

1- SIZE 092 ( or greater )

2- We store these 35 numbers in registers R01 thru R35 like this:

            56   1   2   3   4                  R01  R08  R15  R22  R29
            41   5   1  -3  -1                 R02  R09  R16  R23  R30
            61   4   7   2  -3                 R03  R10  R17  R24  R31
            25   2   3  -1   2      in        R04  R11  R18  R25  R32           respectively.   ( we must store 0  at the place of  F that is in R07 here )
             7    1   1   1   1                 R05  R12  R19   R26  R33
            17   2  -1   1   3                 R06  R13  R20   R27  R34
              0   2   4   3   5                 R07  R17  R21   R28  R35

3- There are 3 inequations of the first type, 1 equation , 2 inequations of the second type and 4 unknowns, therefore:

            3   ENTER^
            1   ENTER^
            2   ENTER^
            4   XEQ "SIMPLEX"  >>>>   F = 76 in the X-register and in R00  ( F01 and F02 are cleared )

            RCL 01 gives  2
            RCL 02   "      7
            RCL 03   "      8
            RCL 04   "      4       thus the maximum of F is 76 and it is achieved for x = 2, y = 7, z = 8, t = 4.

We can also verify that:

            R05 = 0 ; R06 = 52 ; R07 = 0    ( the n slack variables of the three first inequations )
            R08 = R09 = R10 = 0               ( the n'+n" artificial variables which are always 0 when a solution exists )
            R11 = 14 ; R12 = 0                  ( the n" slack variables of the two last inequations )
 
 

NATURAL CUBIC SPLINE INTEGRATION
 

  NCSI  calculates §_x1^xn y dx , given a set of n data points which must be stored as follows:

REGISTERS:

     •  R00 = n = number of points             R2n+1 to R6n-7  are also used

     •  R01 = x₁     •  R03 = x₂      .......  •  R2n-1 = x_n
     •  R02 = y₁     •  R04 = y₂      .......   •  R2n   = y_n

FLAGS:            /
 
 

STACK	INPUT	OUTPUT
X	/	§x1xn  f(x).dx

Example:     Calculate  §₁⁸  f(x).dx   from the following values
 

x	1	2.4	4	5.2	7	8
f(x)	1	4	6	5	4	2

Store these 12 numbers into registers   R01  R03  R05  R07  R09  R11   ( x )
                                                          R02  R04  R06  R08  R10  R12   ( f(x) )    respectively

-SIZE 028   ( SIZE = 6n-8 )

-There are 6 points,     6   STO 00    XEQ "NCSI"    >>>>    §₁⁸  f(x).dx  =  29.99938860
 

TRIDIAGONAL SYSTEMS
 

-"3DLS" is called as a subroutine by "NCSI"
-The coefficients are to be stored in column order, and the system is solved without pivoting

REGISTERS:    R00 & Rbb thru Ree ,  ee =  4n + bbb - 3
FLAGS:            /
 
 

STACK	INPUTS	OUTPUTS
X	(bbb.eee)system	(bbb.eee)solution

   where  (bbb.eee)_system  is the control number of the system    eee =  4n + bbb - 3  , bb > 00 ,  n = number of unknowns ,  n > 1
              (bbb.eee)_solution  --------------------------  solution

Example:

   2.x + 5.y                                     = 2                                                  2   5                           2                 R01   R03                                            R17
   3.x + 7.y + 4.z                            = 4                                                  3   7   4                      4                 R02   R04   R06                                   R18
               y + 3.z + 7.t                    = 7      if you choose  bb = 01               1   3   7                 7                          R05   R07   R09                          R19
                     2.z + 4.t + 6.u           = 1      store these 22 elements                   2   4   6            1          into                    R08   R10   R12                 R20
                              8.t +  u   + 7.v  = 5                                                                 8   1   7       5                                             R11   R13   R15        R21
                                      9.u + 4.v  = 6                                                                      9   4       6                                                       R14   R16        R22

-Then,   1.022  XEQ "3DLS"  >>>>   17.022

  and           R17 = x = -16.47810408            R20 = t = -0.480422463
                  R18 = y =     6.991241632          R21 = u =  0.112313241
                  R19 = z =     1.123905204          R22 = v =  1.247295209

-We have also  R00 = det A = 3882
 

LAGRANGE INTERPOLATION FORMULA
 

-Given a set of n data points:   M₁ ( x₁ , y₁ )  ,   M₂ ( x₂ , y₂ )  ,  ..........  ,  M_n ( x_n , y_n )
  store the 2n numbers   x₁ , y₁ , x₂ , y₂ , ........ , x_n , y_n  into contiguous registers,  starting at any register  Rbb:

     •  Rbb = x₁       •  Rb+2 = x₂      .......     •  Ree-1 = x_n
     •  Rb+1 = y₁     •  Rb+3 = y₂      .......     •  Ree   = y_n                  ( ee = bb + 2n -1 )

  LAGR evaluates  L(x)  from x in X-register and the control number of the array bbb.eee in Y-register.
 

REGISTERS:    Rbb thru Ree ,  ee =  2n + bb - 1
FLAGS:            /
 
 

STACK	INPUTS	OUTPUTS
Y	bbb.eee	bbb.eee
X	x	L(x)
L	/	x

Example:

  Evaluate  y(3) and y(5) for the collocation polynomial that takes the values prescribed below:
 

xi	0	1	2	4	7
yi	3	2	4	6	5

 For instance, you store these 10 numbers ( 0 3 1 2 2 4 4 6 7 5 ) into R01 thru R10, then

 1.010  ENTER^
     3     XEQ "LAGR"  produces:  y(3) = 5.847619048

 RDN    5   R/S           ----------   y(5) = 4.523809520

Note:

-Don't interrupt "LAGR" because synthetic registers P & Q are used
 

INTEGRATION OF THE LAGRANGE POLYNOMIAL
 

  LAGRI  calculates §_x1^x2 L(x) dx , given a set of n data points which must be stored as follows:

     •  R00 = n = number of points

     •  R01 = x₁     •  R03 = x₂      .......  •  R2n-1 = x_n
     •  R02 = y₁     •  R04 = y₂      .......   •  R2n   = y_n

REGISTERS:    R00 thru R2n
FLAGS:            /
 
 

STACK	INPUT	OUTPUT
X	/	§x1xn  L(x).dx

Example:     Calculate  §₁⁸  L(x).dx   from the following values
 

x	1	2.4	4	5.2	7	8
y	1	4	6	5	4	2

Store these 12 numbers into registers   R01  R03  R05  R07  R09  R11   ( x )
                                                          R02  R04  R06  R08  R10  R12    ( y )    respectively
 

-There are 6 points:     6   STO 00      XEQ "LAGRI"     >>>>    §₁⁸  L(x).dx  =   29.61789480

Notes:

-Don't interrupt "BVI" because synthetic registers P & Q are used
-The data registers are altered by this routine.
 

BIVARIATE INTERPOLATION
 

  BVI uses the Lagrange formula  in both axis x & y
  Now we have a grid of  nxm  values:   f(x_i,y_j)   with   i = 1,2,......,n  &  j = 1,2,.....,m   and we want to estimate f(x,y)
  The Data are to be stored as follows:

         •  R00 = n.mmm = n + m/10³

                 x  \  y |       •  Rnn+1 = y₁             •  R2n+2 = y₂       ..............................    •  Rnm+m = y_m
         ---------------------------------------------------------------------------------------------------------------------
        •  R01  = x₁ |   •  Rnn+2 = f(x₁,y₁)    •  R2n+3 = f(x₁,y₂)    ..............................     •  Rnm+m+1 = f(x₁,y_m)
        •  R02  = x₂ |   •  Rnn+3 = f(x₂,y₁)    •  R2n+4 = f(x₂,y₂)    ..............................     •  Rnm+m+2 = f(x₂,y_m)
             ...........    |      .............................................................................................................................................
                            |
         •  Rnn = x_n  |   •  R2n+1  = f(x_n,y₁)   •  R3n+2 = f(x_n,y₂)   ..............................      •  Rnm+m+n  = f(x_n,y_m)
 

REGISTERS:    R00 thru Rn.m+m+n
FLAGS:            /
 
 

STACK	INPUTS	OUTPUTS
Y	y	y
X	x	f(x,y)
L	/	x

Example:

     x \ y |   2    3    4    6                                                                                  |   R04   R08   R12   R16
--------------------------                                                                        --------------------------------
       1   |   4    3    3    5                 are to be stored into                           R01  |  R05   R09   R13   R17            respectively
       2   |   3    1    2    6                                                                           R02  |  R06   R10   R14   R18
       4   |   1    0    4    9                                                                           R03  |  R07   R11   R15   R19

-Here,  n = 3 and  m = 4  so   3.004  STO 00

-Then, to compute  f(3,5)

   5  ENTER^
   3  XEQ "BVI"  >>>>   f(3,5) =  5.833333333

-And to obtain  f(1.6,2.7)

   2.7  ENTER^
   1.6     R/S        >>>>   f(1.6,2.7) =  1.90382

Note:

-Don't interrupt "BVI" because synthetic registers P & Q are used
 

STIRLING NUMBERS OF THE 1st & 2nd KIND
 

-Stirling numbers of the first kind S_n^m are defined by the recurrence relation: ( S_n⁰ = 0 ) ;  S_n^m = S_n-1^m-1 - ( n-1 ) S_n-1^m      ( 1 <= m <= n )
-Stirling numbers of the second kind s_n^m ----------------------------------  ( s_n⁰ = 0 )  ;   s_n^m =  s_n-1^m-1  +  m s_n-1^m            --------------
 

REGISTERS:    SIZE = max ( 2 , k+1 )
FLAG:              F02

   If  CF 02  SNK  computes the Stirling numbers of the 1st kind
   If  SF 02  SNK  ------------------------------------  2nd kind
 
 

STACK	INPUTS	OUTPUTS
Y	n	/
X	k <= n	s(n,k)

where  s(n,k) is the Stirling number of the first kind if flag F02 is clear
                     or ----------------------- second --------  F02 is set.
 

Example:   Calculate  S₁₂⁷   and    s₁₂⁷

   a)  CF 02
        12  ENTER^
          7  XEQ "SNK"   yields       S₁₂⁷ =  -2637558      ( and R01 = -39916800 = S₁₂¹ ; R02 = 120543840 = S₁₂²
                                                                                          .........................................      R07 = -2637558 = S₁₂⁷  )

   b)  SF 02
        12  ENTER^
         7     R/S          produces       s₁₂⁷ = 627396       ( and R01 = 1 =  s₁₂¹ ; R02 = 2047 =  s₁₂² ; ........ ;  R07 = 627396 = s₁₂⁷ )
 
 

CLEBSH-GORDAN COEFFICIENTS & WIGNER 3J-SYMBOL
 

   •  Clebsch-Gordan Coefficients  = <  j₁ j₂  m₁ m₂ |   j₁ j₂ j m  >
                                                    = [ ( 2j +1) D( j₁ j₂ j  ) ]^1/2  [ ( j₁ + m₁ )! ( j₁ - m₁ )! ( j₂ + m₂ )! ( j₂ - m₂ )! ( j + m )! ( j - m )! ]
                                   x Sum_k  (-1)^k / [ k! ( j₁ + j₂ - j - k )! ( j₁ - m₁ - k )! ( j₂ + m₂ - k )! ( j - j₂ + m₁ + k )! ( j - j₁ - m₂ + k )! ]         if   m₁ + m₂ = m

      ( The sum is to be performed for all integers k that produce non-negative integer arguments for the factorials ).

   •  <  j₁ j₂  m₁ m₂ |   j₁ j₂ j m  > = 0            if   m₁ + m₂ # m

   •  Wigner 3j-symbol  = W3J (  j₁ j₂ j  m₁ m₂ -m ) = (-1)^( j₁ - j₂ + m )  ( 2j +1) ^-1/2  <  j₁ j₂  m₁ m₂ |   j₁ j₂ j m  >

  •   Triangular coefficients  D( a b c ) = [ ( a + b - c )! ( a - b + c )! ( -a + b + c )! ] / ( a + b + c + 1 )!
 

>>>  The 6 parameters ( integers or half-integers ) are to be stored in registers R01 thru R06

     •  R01 = j₁     •  R02 = j₂    •  R03 = j
     •  R04 = m₁   •  R05 = m₂   •  R06 = m

REGISTERS:    R00 thru R14
FLAGS:            F10
 
 

STACK	INPUTS	OUTPUTS
Y	/	W3j
X	/	ClG

   where    W3j = W3j (  j₁j₂ j  m₁ m₂ -m )      and  ClG = <  j₁ j₂  m₁ m₂ |   j₁ j₂ j m  >

Examples:

-If   j₁ = 12  ,   j₂ = 24  ,  j  = 31
      m₁ = 1   ,  m₂ = 16  , m = 17

 XEQ "CGW3J"  returns   ClG =  0.206754082
                           X<>Y   W3j = -0.026048566

-With   j₁ = 1  ,   j₂ = 3/2  ,  j  = 5/2
           m₁ = 0 ,  m₂ = 3/2  , m = 3/2

  XEQ "CGW3J"  yields   ClG =  0.632455533
                           X<>Y   W3j = -0.258198890
 

WIGNER 6J-SYMBOL
 

     W6J ( j₁ j₂ j₃ j₄ j₅ j₆ )  =  [ D ( j₁ j₂ j₃ ) D ( j₁ j₅ j₆ ) D ( j₂ j₃ j₄ ) D ( j₃ j₄ j₅ ) ]^1/2
                                       x Sum_{k=k' to k"}  ( -1 )^k ( k + 1 )! / [ ( k - j₁ - j₂ - j₃ )! ( k - j₁ - j₅ - j₆ )! ( k - j₂ - j₃ - j₄ )! ( k - j₃ - j₄ - j₅ )!
                                                                                              ( j₁ + j₂ + j₄ + j₅ - k )! ( j₁ + j₃ + j₄ + j₆ - k )! ( j₂ + j₃ + j₅ + j₆ - k )! ]

       with   k' = max { j₁ + j₂ + j₃ , j₁ + j₅ + j₆ , j₂ + j₃ + j₄ , j₃ + j₄ + j₅ }
       and    k" = min { j₁ + j₂ + j₄ + j₅ , j₁ + j₃ + j₄ + j₆ , j₂ + j₃ + j₅ + j₆ }
 

 Here, we also have 6 numbers ( integers or half-integers ) to store into R01 .... R06

        •  R01 = j₁     •  R02 = j₂    •  R03 = j₃
        •  R04 = j₄     •  R05 = j₅    •  R06 = j₆

REGISTERS:    R00 thru R14
FLAGS:            F10
 
 

STACK	INPUTS	OUTPUTS
X	/	W6j

    with  W6j = W6j ( j₁ j₂ j₃ j₄ j₅ j₆ )

Example:

-If   j₁ = 10  ,   j₂ = 16  ,  j₃  = 21   which are to be stored into      R01   R02   R03
       j₄ = 24  ,   j₅ = 12  ,  j₆  = 14                                                 R04   R05   R06

 XEQ "W6J"  >>>>   W6j  =  0.006251585
 

WIGNER 9J-SYMBOL
 

     W9J ( j₁ j₂ j₃ j₄ j₅ j₆ j₇ j₈ j₉ )  =  Sum_{k=k' to k"}  (-1)^2k (2k+1)  W6J ( j₁ j₉ k j₈ j₄ j₇ )  W6J ( j₈ j₄ k j₆ j₂ j₅ )  W6J ( j₆ j₂ k j₁ j₉ j₃ )
 

    where    k' = max { | j₁ - j₉ | , | j₄ - j₈ | , | j₂ - j₆ | }    and    k" =  min {  j₁ + j₉  ,  j₄ + j₈  ,  j₂ + j₆ }
 

 We now have 9 numbers ( integers or half-integers ) to store into R01 .... R09

       •  R01 = j₁     •  R02 = j₂    •  R03 = j₃
       •  R04 = j₄     •  R05 = j₅    •  R06 = j₆              registers R01 to R09 are saved in registers R15 to R23
       •  R07 = j₇     •  R08 = j₈    •  R09 = j₉              and they are restored at the end

REGISTERS:    R00 thru R27
FLAGS:            F10
 
 

STACK	INPUTS	OUTPUTS
X	/	W9j

  where  W9j = W9J ( j₁ j₂ j₃ j₄ j₅ j₆ j₇ j₈ j₉ )

Examples:

         j₁ = 3  ,   j₂ = 7  ,  j₃  = 5                                                   R01   R02   R03
-If     j₄ = 6  ,   j₅ = 8  ,  j₆  = 9     which are to be stored into      R04   R05   R06
         j₇ = 4  ,   j₈ = 5  ,  j₉  = 7                                                   R07   R08   R09

 XEQ "W9J"  >>>>   W9j  =  0.001032402067

               j₁ = 1/2  ,   j₂ =  1    ,  j₃  = 3/2
-With      j₄ =  1    ,   j₅ = 1/2  ,  j₆  = 1/2
               j₇ = 3/2  ,   j₈ = 1/2  ,  j₉  =  1

 XEQ "W9J"  >>>>   W9j  =  -0.02777777771   ( in fact -1/36 )
 

Notes ( CGW3J W6J W9J ):

-All the j's & m's must be integers or half-integers
-They must satisfy several inequalities so that  FACT is applied to non-negative integers only.
-Otherwise, DATA ERROR will be displayed and the result may be regarded as 0.
 

DOUBLE & TRIPLE INTEGRALS
 

  >>>   DTI   calculates double integrals if flag F03 is clear and triple integrals if F03 is set

  >>>   In both cases, Gauss-Legendre 3-point formula is used and the result is in X-register and in R08
 

REGISTERS:    R00 thru R18 for double integrals, R00 thru R22 for triple integrals.
FLAGS:            F03
 

  •  CF 03 Double Integrals    We want to evaluate:   §_a^b §_u(x)^v(x)  f(x;y) dx dy

      •   R00 = the name of the subroutine which calculates f(x;y)
      •   R01 = a
      •   R02 = b
      •   R03 = n = the number of subintervals into which the intervals (a;b) and (u(x);v(x)) are to be divided.
      •   R04 = the name of the subroutine which calculates u(x).
      •   R05 = the name of the subroutine which calculates v(x).

 >>>  3 subroutines must be keyed into program memory:

      -The first one must take x and y from the X-register and Y-register respectively and calculate f(x;y).
      -The second takes x from the X-register and calculates u(x).
      -The third takes x from the X-register and calculates v(x).
 
 

STACK	INPUT	OUTPUT
X	/	§ab §u(x)v(x) f(x,y).dx.dy

Example:   Evaluate  I = §₁² §_x^{x^2} (1+x⁴ y⁴)^1/2 dx dy.

We must load the 3 needed subroutines into program memory ( any global labels, maximum of 6 characters ),
for instance:

  01  LBL "FF"
  02  *
  03  X^2
  04  X^2
  05  SIGN
  06  LAST X
  07  +
  08  SQRT
  09  RTN
  10  LBL "U"
  11  RTN
  12  LBL "V"
  13  X^2
  14  RTN
 

and then,  "FF" ASTO 00
                   1     STO 01
                   2     STO 02
                 "U"  ASTO 04
                 "V"  ASTO 05      CF 03

  With    n = 1     1 STO 03  XEQ "DTI"  the result is  15.45937082    in the X-register and in R08
             n = 2     2 STO 03      R/S            "      "     "   15.46673275
             n = 4     4 STO 03      R/S            "      "     "   15.46686031
             n = 8     8 STO 03      R/S            "      "     "   15.46686245  ... all the digits are correct!

  •  SF 03 Triple Integrals    We want to evaluate:      §_a^b §_u(x)^v(x) §_w(x;y)^t(x;y)   f (x;y;z) dx dy dz

      •   R00 = the name of the subroutine which calculates f(x;y;z)
      •   R01 = a
      •   R02 = b
      •   R03 = n = the number of subintervals into which the intervals of integration are to be divided.
      •   R04 = the name of the subroutine which calculates u(x).
      •   R05 = -------------------------------------------  v(x).
      •   R06 = ------------------------------------------- w(x;y)
      •   R07 = ------------------------------------------- t(x;y)

>>>  5 subroutines must be keyed into program memory:

      -The first one must take x, y and z from the X-register,the Y-register and the Z-register respectively and calculate f(x;y;z).
      -The second takes x from the X-register and calculates u(x).
      -The third takes x from the X-register and calculates v(x).
      -Another one takes x and y from the X and Y registers respectively and calculates w(x;y)
      -The last one -------------------------------------------------------------------- t(x;y)
 
 

STACK	INPUT	OUTPUT
X	/	§ab §u(x)v(x) §w(x;y)t(x;y)   f (x;y;z) dx dy dz

Example:      Evaluate   I = §₁² §_x^{x^2} §_x+y^xy  xyz  ( x² + y² + z² ) ^-1/2 dx dy dz

  The 5 required subroutines are, for instance:
  ( with global labels of 6 characters maximum )
 

  01  LBL "FF"
  02  ENTER^
  03  X^2
  04  R^
  05  ST* Z
  06  X^2
  07  +
  08  R^
  09  ST* Z
  10  X^2
  11  +
  12  SQRT
  13  /
  14  RTN
  15  LBL "U"
  16  RTN
  17  LBL "V"
  18  X^2
  19  RTN
  20  LBL "W"
  21  +
  22  RTN
  23  LBL "T"
  24  *
  25  RTN

Then we initialize registers R00 thru R07:

    "FF"   ASTO 00
      1        STO 01
      2        STO 02
     "U"    ASTO 04
     "V"    ASTO 05
    "W"    ASTO 06
     "T"    ASTO 07                 SF 03

With  n = 1      1  STO 03    XEQ "DTI"  gives    I₁ =  0.765014888
         n = 2      2  STO 03         R/S            "        I₂ =  0.770640690
         n = 4      4  STO 03         R/S            "        I₄ =  0.770731245
         n = 8      8  STO 03         R/S            "        I₈ =  0.770732669

Note:

-Registers R23 , R24 , .... may be employed to compute
 the functions if they are too much complicated to fit in the stack.
 

-THE EARTH
 

TERRESTRIAL GEODESIC DISTANCE
 

  TGD  uses Andoyer's formula ( second order accuracy ) to compute the geodesic distance between 2 points on the Earth,
            and the forward and backward azimuths:
 

-With   L = ( L₂ - L₁ )/2  ;   F = ( b₁ + b₂ )/2  ;  G = ( b₂ - b₁ )/2
           S = sin² G  cos² L  +  cos² F  sin² L  ;   µ = Arcsin S^1/2  ;  R = ( S.(1-S) )^1/2

 we define  A = ( sin² G  cos² F )/S + ( sin² F  cos² G )/( 1-S )
         and   B = ( sin² G  cos² F )/S - ( sin² F  cos² G )/( 1-S )

-The geodesic distance is then   D = 2a.µ ( 1 + eps )

  where   eps =    (f/2) ( -A - 3B R/µ )
                       + (f²/4) { { - ( µ/R + 6R/µ ).B + [ -15/4 + ( 2S - 1 ) ( µ/R + (15/4) R/µ ) ] A + 4 - ( 2S - 1 ) µ/R } A
                                      - [ (15/2) ( 2S - 1 ).B  R/µ - ( µ/R + 6R/µ ) ] B }

-With the same notations as above, the forward azimuth Az₁ in the direction from P₁ towards P₂ is obtained by the following formulae:    Az₁ = Az_1sph + dAz

      where  dAz = (f/2).( cos² F - sin² G ).[ ( 1+µ/R ).( sin G cos F )/S + ( 1-µ/R ).( sin F cos G )/(1-S) ].( sin 2L )

       and    Tan Az_1sph = cos b₂ sin 2L / ( cos b₁ sin b₂ - sin b₁ cos b₂ cos 2L )     >>>  the R-P function returns the angle in the proper quadrant

  Az_1sph is the forward azimuth in a spherical model of the Earth,
  dAz is a correction that takes the Earth's flattening f into account.

-The back azimuth ( from P₂ to P₁ ) is obtained by replacing L by -L and exchanging b₁ & b₂ in the formulas above.
 

REGISTERS:    R00 thru R09
FLAGS:            /
 
 

STACK	INPUTS	OUTPUTS
T	L1 ( ° ' " )	Az2 ( ° ' " )
Z	b1 ( ° ' " )	Az2 ( ° ' " )
Y	L2 ( ° ' " )	Az1 ( ° ' " )
X	b2 ( ° ' " )	D ( km )

Example:

   -77.0356    ENTER^
    38.55172  ENTER^
     2.20138   ENTER^
    48.50112  XEQ "TGD"  >>>>    D   =   6181.621809 km
                                          RDN    Az₁ =  51°47'36"76 = forward azimuth
                                          RDN    Az₂ = -68°09'59"26 = back azimuth

-The accuracy is usually very good ( errors < 1 meter ) excepy for nearly antipodal points.
 

ELLIPSOIDAL GRAVITY
 

-Let  b = geographic ( geodetic ) latitude and  h = altitude
-The gravity g above an ellipsoid of revolution may be computed by the rigorous formulae:

             g  =  ( p²+ q² )^1/2

 where   p =  { - [ GM + ( omega )² a² E ( q' / 2q₀ ) ( sin² F - 1/3 ) ] / ( u² + E² ) + ( omega )² u cos² F } / w
  and     q =   [ a² ( q / q₀ ) ( u² + E² ) ^-1/2 - ( u² + E² )^1/2 ] ( omega )² ( sin F cos F ) / w

            u  =  { [ ( x² + y² + z² - E² ) / 2 ] [ 1 + { 1 + 4 E² z² / ( x² + y² + z² - E² )² }^1/2 ] }^1/2
  with    E =   a ( 2 f - f ² )^1/2  ,   sin F = z / u   ,   cos F = ( x² + y² )^1/2 / ( u² + E² )^1/2
            w =  [ ( u² + E² sin² F )^1/2 / ( u² + E² )^1/2 ]^1/2

            q  = (1/2) [ ( 1 + 3 u² / E² ) Atan ( E / u )  - 3 E / u ]                                                            ( x² + y² )^1/2  = ( N + h ) cos b
  and    q₀ = (1/2) [ ( 1 + 3 b' ² / E² ) Atan ( E / b' ) - 3 E / b' ]         ( b' = polar radius )                          z            = [ ( b'/a )² N + h ] sin b
            q' =   3 ( 1 + u² / E² ) [ 1 - ( u / E ) Atan ( E / u ) ] - 1                                                                 N          = a ( 1 - e² sin² b ) ^-1/2
 

-If, however, we use ATAN to calculate  q , q₀ and q' , the results will not be accurate because of cancellation of leading digits.
-So, the following program employs ascending series to compute:

       Atan t + 3 [ ( Atan t ) / t - 1 ] / t      =  4 t³  Sum_{k=0,1,2,.....}  ( k + 1 ) ( - t² )^k / [ ( 2 k + 3 ) ( 2 k + 5 ) ]
    3 ( 1 + 1 / t² ) [ 1 - ( Atan t ) / t ] - 1  =  6 t²  Sum_{k=0,1,2,.....}  ( - t² )^k / [ ( 2 k + 3 ) ( 2 k + 5 ) ]

-Since the argument t never exceeds 0.0821 for the Earth, only a few terms must be calculated.
 

REGISTERS:    R00 thru R13
FLAGS:            F10
 

  GRV  takes the latitude ( in deg.mnss ) in Y-register and the altitude ( in meters ) in X-register,
            it returns the gravity ( in m/s² ) in X
            the distance from the center of the Earth ( in meters ) in Y
            and the geocentric latitude ( in deg )  in Z
 
 

STACK	INPUTS	OUTPUTS
Z	/	B ( deg )
Y	b ( ° ' " )	R ( m )
X	h ( m )	g ( m/s2 )

   where   R = distance from the center of the Earth
    and      B = geocentric latitude

Example1:     b = 38°55'17"2     h = 23456 m

  38.55172   ENTER^
    23456      XEQ "GRV"  >>>>   g = 9.728750374 m/s²
                                         RDN    R = 6393195.112 m
                                         RDN    B = 38°73415897
 

Example2:     b = 38°55'17"2     h = 12345678 m

  38.55172   ENTER^
 12345678       R/S      >>>>   g = 1.078713288 m/s²
                                    RDN   R = 18715394.62 m
                                    RDN   B = 38°85746766

Note:

-The formula is rigorous for an ellipsoid of revolution
-For the Earth, this is only a ( good ) approximation ...
 

EULER's GAMMA CONSTANT
 

-The M-Code routine GEU  simply "recalls"  Euler's Constant = 0.5772156649  in X-register
-Register L is unchanged.
 
 

STACK	INPUTS	OUTPUTS
T	t	z
Z	z	y
Y	y	x
X	x	0.5772156649

 

-SPEC FCNS
 

AIRY FUNCTIONS
 

    Ai(x) = [ 3 ^-2/3 / Gam(2/3) ] ₀F₁( ; 2/3 ; x³/9 ) - x [ 3 ^-1/3 / Gam(1/3) ] ₀F₁( ; 4/3 ; x³/9 )
    Bi(x) = [ 3 ^-1/6 / Gam(2/3) ] ₀F₁( ; 2/3 ; x³/9 ) + x [ 3 ^1/6 / Gam(1/3) ] ₀F₁( ; 4/3 ; x³/9 )
 

REGISTERS:    R00 thru R04
FLAGS:            /
 
 

STACK	INPUTS	OUTPUTS
Y	/	Bi(x)
X	x	Ai(x)

Example:

  0.4  XEQ "AIRY"  >>>>   Ai(0.4) = 0.254742355
                               X<>Y   Bi(0.4) = 0.801773001
 

ERROR FUNCTION
 

-This function is defined by:        erf x =  (2/pi^1/2)  §₀^x  e^{-t^2} dt

      Ascending series                 erf x = (2x/Pi^1/2) exp(-x²)  ₁F₁( 1 , 3/2 , x² )

    Continued Fraction:      2.e^{x^2}  §_x^infinity e^{-t^2} dt = 1/(x+0.5/(x+1/(x+1.5/(x+2/(x+....)))))       the second formula is used if x > Ln 6
 

REGISTERS:    R00 thru R02
FLAGS:            /
 
 

STACK	INPUTS	OUTPUTS
Y	/	1-erf x
X	x >= 0	erf x

Examples:         0.9    XEQ ERF  >>>>   erf(0.9) = 0.796908213      X<>Y     1-erf(0.9) = 0.203091787
                          2.7    XEQ ERF  >>>>   erf(2.7) = 0.999865667      X<>Y     1-erf(2.7) = 0.0001343327399
 
 

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) = ( PI x³/6 ) ₁F₂( 3/4 ; 3/2 , 7/4 ; -PI² x⁴/16 )

                        C(x) = x ₁F₂( 1/4 ; 1/2 , 5/4 ; -PI² x⁴/16 )

-When x > sqrt(Pi) , a complex continued fraction is employed:

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

REGISTERS:    R00 thru R04
FLAGS:            /
 
 

STACK	INPUTS	OUTPUTS
Y	/	S(x)
X	x	C(x)

Examples:       1.5  XEQ "CSX"  >>>>   C(1.5) = 0.445261176      X<>Y     S(1.5) = 0.697504960
                          4    XEQ "CSX"  >>>>     C(4)  = 0.498426033      X<>Y       S(4)  = 0.420515754
 
 

SINE INTEGRAL
 

   Si(x)  = §₀^x sin(t)/t dt = x ₁F₂( 1/2 ; 3/2 , 3/2 ; -x²/4 )
 

REGISTERS:    R00 thru R03
FLAGS:            /
 
 

STACK	INPUT	OUTPUT
X	x	Si(x)

     1.4   XEQ "SI"     >>>>  Si(1.4) =  1.256226733
 

HYPERBOLIC SINE INTEGRAL
 

  Shi(x) = §₀^x sinh(t)/t dt = x ₁F₂( 1/2 ; 3/2 , 3/2 ; x²/4 )
 

REGISTERS:    R00 thru R03
FLAGS:            /
 
 

STACK	INPUT	OUTPUT
X	x	Shi(x)

     1.4   XEQ "SHI"  >>>> Shi(1.4) =  1.561713390
 

COSINE INTEGRAL
 

  Ci(x) = C + Ln(x) + §₀^x (cos(t)-1)/t dt = C + Ln(x) - (x²/4) ₂F₃( 1 , 1 ; 2 , 2 , 3/2 ; -x²/4 )               where  C = 0.5772156649 is Euler's constant.
 

REGISTERS:    R00 thru R05
FLAGS:            /
 
 

STACK	INPUT	OUTPUT
X	x	Ci(x)

     1.4   XEQ "CI"    >>>>  Ci(1.4) =  0.462006585
 

HYPERBOLIC COSINE INTEGRAL
 

  Chi(x)= C + Ln(x) + §₀^x (cosh(t)-1)/t dt = C + Ln(x) + (x²/4) ₂F₃( 1 , 1 ; 2 , 2 , 3/2 ; x²/4 )               where  C = 0.5772156649 is Euler's constant.
 

REGISTERS:    R00 thru R05
FLAGS:            /
 
 

STACK	INPUT	OUTPUT
X	x	Chi(x)

     1.4   XEQ "CHI" >>>> Chi(1.4) =  1.445494076
 

COSINE & SINE INTEGRALS  for  X > 3
 

-For x > 3 , a complex continued fraction is preferable:

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

REGISTERS:    R00-R01
FLAGS:           /
 
 

STACK	INPUTS	OUTPUTS
Z	/	Si(x)-pi/2
Y	/	Si(x)
X	x >= 3	Ci(x)

  3  XEQ "CSIX>3"  >>>>        Ci(3)  = 0.119629786
                                 RDN         Si(3)  = 1.848652528
                                 RDN  Si(3)-pi/2  = 0.277856201
 

EXPONENTIAL INTEGRAL Ei(x)
 

   Ei(x) = §_-infinity^x  e^t/t dt  = C + Ln(x) + x ₂F₂( 1 , 1 ; 2 , 2 ; x )                 where  C = 0.5772156649 is Euler's constant.
 

REGISTERS:    R00 thru R05
FLAGS:            /
 
 

STACK	INPUT	OUTPUT
X	x	Ei(x)

     1.4   XEQ "EI"     >>>>    Ei(1.4) =  3.007207464
 

EXPONENTIAL INTEGRAL E_n(x)
 

-We have  E_n(x) = §₁^+infinity  t ^-n e ^-x.t dt        ( x > 0 ;  n = a positive integer )
-This function is computed by a series expansion if x <= 1.5

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

 and by the following continued fraction if  x > 1.5

     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(x) = x^n-1 Gam(1-n) - [1/(1-n)] ₁F₁( 1 , 2-n ; x ) e ^-x

-This last relation is employed if n is not an integer and shouldn't be used for large x-values.
 

REGISTERS:    R00 thru R04
FLAGS:            /
 
 

STACK	INPUT	OUTPUT
Y	n	/
X	x	En(x)

Examples:

       0    ENTER^
     1.4  XEQ "ENX"  >>>>  E₀(1.4) = 0.176140689

       2    ENTER^
     1.6  XEQ "ENX"  >>>>  E₂(1.6) = 0.063803184

    -2.1   ENTER^
     3.4   XEQ "ENX"  >>>>  E_-2.1(3.4) = 0.017887050
 

RIEMANN ZETA FUNCTION
 

-The Riemann Zeta function is defined by the series:       Zeta(x) = 1 + 1/2^x + 1/3^x + ........ + 1/n^x + ......     ( for x > 1)
  but the following program uses the formula:          Zeta(x) = ( 1 - 1/2^x + 1/3^x - 1/4^x + ...... ) / ( 1 - 2^1-x )   where roundoff errors are smaller.

-If  x < 1, we have:   Zeta(x) = GAM(1-x)  2^x  PI^x-1  Sin(PI.x/2)  Zeta(1-x)
 

REGISTERS:    R00 thru R03
FLAGS:            /
 
 

STACK	INPUT	OUTPUT
X	x	Zeta(x)

Examples:

  FIX 7       3    XEQ "ZETA"    returns    Zeta(3) = 1.2020569
  FIX 9       3    XEQ "ZETA"    -------    Zeta(3) = 1.202056903

  FIX 9   -7.49  XEQ "ZETA"   -------   Zeta(-7.49) = 0.003312040169

Notes:

-Execution time and Zeta(x) tend to infinity as x tends to 1.
-The accuracy is determined by the display setting.
 

MODIFIED BESSEL FUNCTION 1st KIND

                             ~
    I_n(x) = ( x/2 )ⁿ ₀F₁( ; n+1 ; x²/4 )
 

REGISTERS:    R00 thru R03
FLAGS:            F01
 
 

STACK	INPUTS	OUTPUTS
Y	n	/
X	x	In(x)

   0.7   ENTER^
   1.9   XEQ "INX"   yields   I_0.7(1.9) = 1.727630603

-The results are accurate for large arguments too.
 

BESSEL FUNCTION 1st KIND

                            ~
   J_n(x) = ( x/2 )ⁿ ₀F₁( ; n+1 ; - x²/4 )
 

REGISTERS:    R00 thru R03
FLAGS:            F01
 
 

STACK	INPUTS	OUTPUTS
Y	n	/
X	x	Jn(x)

   0.7   ENTER^
   1.9   XEQ "JNX"   >>>>   J_0.7(1.9) = 0.584978103

-The results are not accurate for large arguments: use "AEJYNX" in this case.
 

BESSEL FUNCTION 1st KIND - NON-NEGATIVE INTEGER ORDER
 

 "JNX1" uses Miller algorithm to produce accurate results for large x-values.

-It uses the relations:   J₀(x) + 2 J₂(x) + 2 J₄(x) + 2 J₆(x) + ...... = 1  and   J_n-1(x) + J_n+1(x) = (2n/x) J_n(x)
 

REGISTERS:    R00 thru R05
FLAGS:            /
 
 

STACK	INPUTS	OUTPUTS
Y	n	/
X	x	Jn(x)

     3    ENTER^
   100  XEQ "JNX1"  >>>>    J₃(100)  =  7.628420178 10 ^-2
 

BESSEL FUNCTION 2nd KIND
 

  Y_n(x) = ( J_n(x) cos(n(pi)) - J_-n(x) ) / sin(n(pi))    if n is not an integer and otherwise ( n non-negative ):

  Y_n(x) = -(1/pi) (x/2)^-n SUM_{k=0,1,....,n-1}  (n-k-1)!/(k!) (x²/4)^k + (2/pi) ln(x/2) J_n(x) - (1/pi) (x/2)ⁿ SUM_k=0,1,.....  ( psi(k+1) + psi(n+k+1) ) (-x²/4)^k / (k!(n+k)!)
 

REGISTERS:    R00 thru R05
FLAGS:            F01
 
 

STACK	INPUTS	OUTPUTS
Y	n	/
X	x	Yn(x)

Examples

    1.4  ENTER^
     3   XEQ "YNX"   >>>>   Y_1.4(3) = 0.137821837

       2   ENTER^
       3   XEQ "YNX"   >>>>   Y₂(3) = -0.160400393

-For large arguments, "AEJYNX" is preferable...
 

MODIFIED BESSEL FUNCTION 2nd KIND
 

      K_n(x) = (pi/2) ( I_-n(x) - I_n(x) ) / sin(n(pi))      if n is not an integer and otherwise ( n non-negative ):

K_n(x) =  (1/2) (x/2)^-n SUM_{k=0,1,....,n-1}  (n-k-1)!/(k!) (-x²/4)^k - (-1)ⁿ ln(x/2) I_n(x) + (1/2) (-1)ⁿ (x/2)ⁿ SUM_k=0,1,.....  ( psi(k+1) + psi(n+k+1) ) (x²/4)^k / (k!(n+k)!)
 

REGISTERS:    R00 thru R05
FLAGS:            F01
 
 

STACK	INPUTS	OUTPUTS
Y	n	/
X	x	Kn(x)

Examples

      1.4  ENTER^
       3    XEQ "KNX"    >>>>   K_1.4(3) = 0.046088044

       2   ENTER^
       3   XEQ "KNX"    >>>>    K₂(3) =  0.061510458

-For large arguments, "AEKNX" is preferable...
 

PARABOLIC CYLINDER FUNCTION
 

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

REGISTERS:    R00 thru R05
FLAGS:            /
 
 

STACK	INPUTS	OUTPUTS
Y	n	/
X	x	Dn(x)

   0.4   ENTER^
   1.8   XEQ "DNX"  >>>>   D_0.4(1.8) = 0.579579485

-For large x-values, "AEDNX" is preferable...
 

KELVIN FUNCTION 1st KIND
 

Formulae:

   ber_n(x) = [ cos(135°n) / Gam(n+1) ] (x/2)ⁿ S  -  [ sin(135°n) / Gam(n+2) ] (x/2)ⁿ⁺² T

   bei_n(x) = [ cos(135°n) / Gam(n+2) ] (x/2)ⁿ⁺² S  +  [ sin(135°n) / Gam(n+1) ] (x/2)ⁿ T
 

 where      S = ₀F₃( ; (n+1)/2 , 1+n/2 ,1/2 ; - x⁴/256 )       T = ₀F₃( ; 1+n/2 , (n+3)/2 , 3/2 ; - x⁴/256 )
 

REGISTERS:    R00 thru R07
FLAGS:            /
 
 

STACK	INPUTS	OUTPUTS
Y	n	bein(x)
X	x	bern(x)

  where n is not a negative integer.

Example:

   2   SQRT
  PI   XEQ "KLV"  >>>>    ber_sqrt(2)(PI) = -0.674095953
                            X<>Y    bei_sqrt(2)(PI) = -1.597357212

-For large arguments, "AEKLV" is preferable.
 

KELVIN FUNCTION 2nd KIND
 

-If n is not an integer, we have:

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

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

-Otherwise,

     ker_n(x) = (1/2) (x/2)^-n SUM_{k=0,1,....,n-1} cos ((3n/4+k/2)180°) (n-k-1)!/(k!) (x²/4)^k - ln(x/2) ber_n(x) + (PI/4) bei_n(x)
                 + (1/2) (x/2)ⁿ SUM_k=0,1,.....   cos ((3n/4+k/2)180°)  ( psi(k+1) + psi(n+k+1) ) (x²/4)^k / (k!(n+k)!)

     kei_n(x) = - (1/2) (x/2)^-n SUM_{k=0,1,....,n-1} sin ((3n/4+k/2)180°) (n-k-1)!/(k!) (x²/4)^k - ln(x/2) bei_n(x) - (PI/4) ber_n(x)
                   + (1/2) (x/2)ⁿ SUM_k=0,1,.....   sin ((3n/4+k/2)180°)  ( psi(k+1) + psi(n+k+1) ) (x²/4)^k / (k!(n+k)!)

  where  psi = digamma function
 

REGISTERS:    R00 thru R10
FLAGS:            /
 
 

STACK	INPUTS	OUTPUTS
Y	n	kein(x)
X	x	kern(x)

  where n is not a negative integer.

Examples:

   2   SQRT
  PI   XEQ "KLV2"  >>>>    ker_sqrt(2)(PI) =  0.025901894
                               X<>Y    kei_sqrt(2)(PI) =  0.089242866

   3    ENTER^
  PI    XEQ "KLV2"    >>>>    ker₃(PI) =  -0.045754846
                                  X<>Y    kei₃(PI) =  -0.196928536

-For large arguments, "AEKLV" is preferable.
 

STRUVE FUNCTION

                                 ~
    H_n(x) =  (x/2)ⁿ⁺¹ ₁F₂( 1 ; 3/2 , n+3/2 ; - x²/4 )
 

REGISTERS:    R00 thru R04
FLAGS:            F01
 
 

STACK	INPUTS	OUTPUTS
Y	n	/
X	x	Hn(x)

   1.2   ENTER^
   3.4   XEQ "HNX"  >>>>   H_1.2(3.4) = 1.113372658

-Not to be confused with Hermite functions.
-Use "AEHNX" for large x-values.
 

MODIFIED STRUVE FUNCTION

                                 ~
   L_n(x) =    (x/2)ⁿ⁺¹₁F₂( 1 ; 3/2 , n+3/2 ; x²/4 )
 

REGISTERS:    R00 thru R04
FLAGS:            F01
 
 

STACK	INPUTS	OUTPUTS
Y	n	/
X	x	Ln(x)

     1.2   ENTER^
     3.4   XEQ "LNX"  >>>>   4.649129471

-The output is accurate for large x-values.
 

BETA FUNCTION
 

     B(x,y) = B(y,x) = GAM(x).GAM(y)/GAM(x+y)
 

REGISTERS:    /
FLAGS:            /
 
 

STACK	INPUTS	OUTPUTS
Y	y	/
X	x	B(x,y)

    1  E^X   PI   XEQ "BETA"  >>>>  B(e,Pi) = 0.037890299
 

BETA FUNCTION - INTEGER ARGUMENTS
 

    B(m,n) = (m-1)!(n-1)!/(m+n-1)!
 

REGISTERS:    /
FLAGS:            /
 
 

STACK	INPUTS	OUTPUTS
Y	m	/
X	n	B(m,n)

  100  ENTER^
  200  XEQ "BMN"  >>>>   B(100,200) = 3.607285453 10 ^-84
 

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)
 

REGISTERS:    R00 thru R04
FLAGS:            /
 
 

STACK	INPUTS	OUTPUTS
Z	a	/
Y	b	/
X	x	Bx(a,b)

    PI   ENTER^
     1   E^X
   0.7  XEQ "IBF"  >>>>  B_0.7(Pi,e) =  0.029623046
 

INCOMPLETE GAMMA FUNCTION
 

-The incomplete gamma function is defind by    igam(a,x) =  §₀^x  e ^-t t^a-1 dt
  but this program uses the following relation:

     igam(a,x) =  (x^a/a).exp(-x) ₁F₁(1,a+1,x)    where    ₁F₁ = M = Kummer's function.
 

REGISTERS:    R00 thru R03
FLAGS:            /
 
 

STACK	INPUTS	OUTPUTS
Y	a	/
X	x	igam(a,x)

     1.2  ENTER^
     1.7  XEQ "IGF"      >>>>   igam( 1.2 , 1.7 ) = 0.697290898
 

WEBER & ANGER FUNCTIONS
 

-Anger's function is defined by     J_n(x) = (1/PI) §₀^pi  cos ( n.t - x.sin t ) dt
-Weber's function is defined by    E_n(x) = (1/PI) §₀^pi  sin ( n.t - x.sin t ) dt

-The following programs use the formulae:
                                               ~                                                                         ~
   J_n(x) = + (x/2) sin ( 90°n )  ₁F₂( 1 ; (3-n)/2 , (3+n)/2 ; -x²/4 )  + cos ( 90°n ) ₁F₂( 1 ; (2-n)/2 , (2+n)/2 ; -x²/4 )
                                                ~                                                                       ~
   E_n(x) = - (x/2) cos ( 90°n )  ₁F₂( 1 ; (3-n)/2 , (3+n)/2 ; -x²/4 ) + sin ( 90°n ) ₁F₂( 1 ; (2-n)/2 , (2+n)/2 ; -x²/4 )
 

REGISTERS:    R00 thru R06
FLAGS:            /
 
 

STACK	INPUTS	OUTPUTS
Y	n	Jn(x)
X	x	En(x)

   2   SQRT
  PI  XEQ "WEBAN"  >>>>    E_sqrt(2)(PI) = -0.315594385
                                  X<>Y    J_sqrt(2)(PI) =  0.366086559
 

DEBYE FUNCTION
 

   dby(x;n) = §_x^+infinity  tⁿ/(e^t-1).dt

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

REGISTERS:    R00 thru R03
FLAGS:            /
 
 

STACK	INPUTS	OUTPUT
Y	n	/
X	x	§x+infinity  tn/(et-1).dt

     n = positive integer ,  x > 0

Example:

    3   ENTER^
  0.7  XEQ "DBY"  >>>>  dby( 0.7 ; 3 ) = 6.406833597
 

ASSOCIATED LEGENDRE FUNCTIONS 1st KIND

                                            ~
    P_n^m(x) = |(x+1)/(x-1)|^m/2₂F₁(-n ; n+1 ; 1-m ; (1-x)/2 )        (  x # 1 )
 

REGISTERS:    R00 thru R05
FLAGS:            /
 
 

STACK	INPUTS	OUTPUTS
Z	m	/
Y	n	/
X	x	Pnm(x)

Examples:

        0.4   ENTER^
        1.3   ENTER^
        0.7   XEQ "ALF"  >>>>   P_1.3^0.4(0.7) =  0.274932821

   -0.6   ENTER^
    1.7   ENTER^
    4.8   XEQ "ALF"  >>>>>  P_1.7^-0.6(4.8) = 10.67810281
 

ASSOCIATED LEGENDRE FUNCTIONS 2nd KIND
 

Formulae:

   Q_n^m(x) = 2^m pi^1/2  (1-x²)^-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 ; x² )
                                                                     + 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 ; x² )  ]

   Q_n^m(x) = e^i.mpi 2^-n-1 pi^1/2 x^-n-m-1 (x²-1)^m/2 ₂F₁(1+n/2+m/2 ; 1/2+n/2+m/2 ; n+3/2 ; 1/x² )  Gamma(1+n+m) / Gamma(n+3/2)    if  x > 1

    ( the result may be a complex number )
 

REGISTERS:    R00 thru R09
FLAGS:            /
 
 

STACK	INPUTS	OUTPUTS
Z	m	/
Y	n	/
X	x	Qnm(x)

     0.4  ENTER^
     1.3   ENTER^
     0.7   XEQ "ALF2"   >>>>     Q_1.3^0.4(0.7) =  -1.317935681

>>>  If  x > 1 , the HP-41 displays "X+I.Y" and the result is usually a complex number ( if Y # 0 )

    0.7   ENTER^
    1.2   ENTER^
    1.9   XEQ "ALF2"  >>>>    -0.081513570
                                  X<>Y     0.112193805     thus,     Q_1.2^0.7(1.9) = -0.081513570 + 0.112193805 i
 

F-DISTRIBUTION FUNCTION
 

  FDF   calculates    Q( x | n₁ , n₂ ) =  (n₁/n₂)^(n1)/2  Gamma[(n₁+n₂)/2] / [ Gamma(n₁/2) Gamma(n₂/2) ] §_x^infinity  t ^-1+(n1)/2 ( 1 + t.n₁/n₂ ) ^-(n1+n2)/2 dt

  via the incomplete beta function:   Q( x | n₁ , n₂ ) = I_y( n₂/2 , n₁/2 )   with  y = n₂/(n₂+n₁.x)
 

REGISTERS:    R00 thru R05
FLAG:              F09
 
 

STACK	INPUTS	OUTPUTS
Z	n1	/
Y	n2	/
X	x >= 0	Q( x | n1 , n2 )

Examples:

     7    ENTER^
     6    ENTER^
   4.3   XEQ "FDF"  >>>>    Q ( 4.3 |  7 , 6 ) = 0.047640800

     7    ENTER^
     9    ENTER^
   0.4   XEQ "FDF"  >>>>    Q ( 0.4 | 7 , 9 ) =  0.879873501
 

CHI2-PROBABILITY FUNCTION
 

   PX2N computes  P( X² | m ) = 2 ^-m/2 [ 1/Gamma(m/2) ] §₀^{X^2}   t ^-1+m/2 e ^-t/2 dt

-It uses the series exansion:    P( X² | m ) = (X²/2)^m/2 e ^{-(X^2)/2} [ 1/Gamma(m/2) ]  [ 1 + Sum_k=1,2,...  X^2k/[ (m+2)(m+4) ..... (m+2k) ]
 

REGISTERS:    R00 thru R02
FLAGS:            /
 
 

STACK	INPUTS	OUTPUTS
Y	m	/
X	X2	P( X2 | m )

  where m is a positive integer.

Example:

    7   ENTER^
   12  XEQ "PX2N"  >>>>   P (12 | 7 ) = 0.899441132
 

STUDENT's t-DISTRIBUTION FUNCTION
 

- "STD"  computes  A( t | n ) = n ^-1/2 / Beta(0.5 ; n/2 )  §_-t^+t  ( 1 + x²/n ) ^-(n+1)/2  dx

-This routine also computes   n ^-1/2 / Beta(0.5 ; n/2 )  §_t^+infinity  ( 1 + x²/n ) ^-(n+1)/2  dx = [ 1 - A ( t | n ) ] / 2 = UTPT(t,n)

Formula:   A( t | n ) = 1 - I_y( n/2 , 1/2 )   with   y = n/(n+t²)
 

REGISTERS:    R00 thru R05
FLAG:              F09
 
 

STACK	INPUTS	OUTPUTS
Y	n	UTPT(t,n)
X	t >= 0	A( t | n )

Examples:

   4   ENTER^
   1   XEQ "STD"   >>>>            A( 1 | 4 ) = 0.626099033
                              RDN     UTPT( 1 , 4 ) = 0.186950483

  4   ENTER^
 10  XEQ "STD"    >>>>            A( 10 | 4 ) = 0.999437997
                              RDN     UTPT( 10 , 4 ) = 0.0002810018113
 

REGULAR COULOMB WAVE FUNCTION
 

Formulae:         F_L(n,r) = C_L(n) r ^L+1 Sum _k>L A_k^L (n) r ^k-L-1

    with   C_L(n) = (1/Gam(2L+2)) 2^L e ^-pi.n/2 | Gam(L+1+i.n) |
    and    A_L+1^L = 1 ;  A_L+2^L = n/(L+1)  ;  (k+L)(k-L-1) A_k^L = 2n A_k-1^L - A_k-2^L   ( k > L+2 )

-We also use   | Gam( 1+i y ) |² = (Pi.y) / Sinh (Pi y)   and

    | Gam( 1+L+i y ) |² = [ L² + y² ] [ (L-1)² + y² ] .................. [ 1 + y² ]  (Pi.y) / Sinh (Pi y)
 

REGISTERS:    R00 thru R06
FLAGS:            F10
 
 

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 non-negative.

Example:

     2   ENTER^
   0.7  ENTER^
   1.8  XEQ "RCWF"  >>>>   F₂( 0.7 , 1.8 ) = 0.141767746
 

TORONTO FUNCTION
 

Toronto function is defined by

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

REGISTERS:    R00 thru R04
FLAGS:            /
 
 

STACK	INPUTS	OUTPUTS
Z	m	/
Y	n	/
X	r	T(m,n,r)

Example:

   2    SQRT
   3    SQRT
  PI   XEQ "TMNR"  >>>>  T( sqrt(2),sqrt(3),PI ) = 0.963524225
 

SPHEROIDAL EIGENVALUES
 

  "LMNC2"   finds 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)) ]
 

REGISTERS:    R00 thru R09
FLAGS:            F01
 
 

STACK	INPUTS	OUTPUTS
Z	m	/
Y	n	/
X	c2	Lmn(c2)

Example1:      m = 4 ,  n = 11 ,  c² = -1    ( note that  c²  may be positive or negative )

    4   ENTER^
   11  ENTER^
   -1  XEQ "LMNC2"   >>>>   L_mn = 131.5600809

Example2:      m = 0.2  ,  n = 0.6 ,  c² = 1.7        ( actually, m & n are not necessarily integers )

   0.2   ENTER^
   0.6   ENTER^
   1.7   XEQ "LMNC2"    >>>>   L_mn = 2.246866650

Notes:

-Usually, m & n are non-negative integers and n = m , m+ 1 , m + 2 , ........
  c² = 0 corresponds to Legendre functions and in that case, L = n ( n + 1 )

-However, L may also be defined mathematically for more general arguments.
-For some values, the algorithm diverges and no solution is returned.
-Avoid n = -1/2

-Though R00 = 12 seems enough, the number of terms in the continued fraction ( half of this value ) is not known in advance,
  so you could press  XEQ 01 after storing - say 14 - in R00
-If the last decimal is ( almost ) the same, the result should be correct.

-The 4 coefficients m , n , c² , L  are automatically stored in R01 , R02 , R03 , R04 at the end, as required by the following program.
 

ANGULAR SPHEROIDAL WAVE FUNCTION of the 1st KIND
 

  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
 

REGISTERS:    R00 thru R10
FLAGS:           F10
 
 

STACK	INPUT	OUTPUT
X	x	Smn(x)

  Provided that m , n , c² , L are stored in R01 thru R04 - which is already done if you have computed L with "LMNC2"

Example:     m = 0.2   n = 0.6   c² = 1.7    "LMNC2"  gives   2.246866650       -Store these 4 numbers into  R01 thru R04  and

    0.7   XEQ "SWF"  >>>>     S_mn(0.7) = 0.682645660

Notes:

-Here,  -1 <= x <= +1
-"SWF" will work even if L does not satisfy the equation that "LMNC2" solves.
-In this case however, the solution is not regular at  x = +/- 1
 

THE HYPERGEOMETRIC FUNCTION
 

  HGF calculates  ₂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)

-The relation        ₂F₁(a,b,c,x) =  ( 1 - x ) ^-a  ₂F₁(a,c-b,c,-x/(1-x))      is used if  x < 0
 

REGISTERS:    R00 thru R03

  •   R01 = a
  •   R02 = b              registers R01 R02 R03 are to be initialized before executing "HGF"
  •   R03 = c

FLAGS:            /
 
 

STACK	INPUTS	OUTPUTS
X	x < 1*	2F1(a,b,c,x)

 * unless F is a polynomial.

Examples:

  •   1.2  STO 01    2.3  STO 02    3.7  STO 03

     0.4  XEQ "HGF"  >>>>   1.435242953
     -3   XEQ "HGF"  >>>>   0.309850661

Note:

-Press any key to stop a possible infinite loop if the series is divergent ( or converges too slowly ).
 

THE REGULARIZED HYPERGEOMETRIC FUNCTION

                               ~
  RHGF calculates  ₂F₁(a,b,c,x)
 

REGISTERS:    R00 thru R03

  •   R01 = a
  •   R02 = b              registers R01 R02 R03 are to be initialized before executing "RHGF"
  •   R03 = c

FLAGS:            /
 
 

STACK	INPUTS	OUTPUTS
X	x < 1*	2F1~(a,b,c,x)

              ~
 * unless F is a polynomial.

Example:

  •   2  STO 01    3  STO 02    -7  STO 03

     0.4  XEQ "RHGF"  >>>>   5353330.290
     -3   XEQ "RHGF"  >>>>   2128.650875

Note:

-Press any key to stop a possible infinite loop if the series is divergent ( or converges too slowly ).
 

-ASYMPTOT !
 

-The following routines - except "AEKLV" - use asymptotic expansions via the M-Code function HGF+
-Since these series are actually divergent, press any key to stop the infinite loop if x is not large enough.
-In this case, the results will be probably meaningless...

-So, there is sometimes a transition region where the asymptotic expansions cannot be used,
 whereas the ascending series that are employed in the main programs return low accurate results.
 

AIRY FUNCTIONS
 

With x > 0

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

With  x > 0    Ai(-x) ~ ( x ^-1/4 / sqrt(PI) ) [ S sin ( z + 45° ) - T cos ( z + 45° ) ]
                      Bi(-x) ~ ( x ^-1/4 / sqrt(PI) ) [ S cos ( z + 45° ) + T sin ( z + 45° ) ]

    where   z = (2/3) x^3/2    ;     S = Sum_k=0,1,2,.... (-1)^k c_2k z ^-2k    ,    T = Sum_k=0,1,2,.... (-1)^k c_2k+1 z ^-2k-1

       and       c₀ = 1 ,  c_k = c_k-1 ( 6k -5 )( 6k - 1 ) / (72 k )
 

REGISTERS:    R00 thru R07
FLAGS:            /
 
 

STACK	INPUTS	OUTPUTS
Y	/	Bi(x)
X	x	Ai(x)

Examples:

  6.4   XEQ "AEAIRY"  >>>>   Ai(6.4) = 0.000003617762307
                                       RDN   Bi(6.4) = 17400.13559

-With x = 6.3 ,  the series diverge too soon:  press any key to stop the infinite loop.

 -7.4   XEQ "AEAIRY"  >>>>   Ai(-7.4) =   0.341323752
                                      RDN    Bi(-7.4) =  -0.021596519

-With x = -7.3 ,  the series diverge too soon.
 

BESSEL FUNCTION 1st KIND & 2nd KIND
 

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 =  ₄F₁( (1-2n)/4 , (1+2n)/4 ,  (3-2n)/4 , (3+2n)/4 ; 1/2 ; -1/x² )
     and   Q =  [(1-4n²)/(8x)] ₄F₁( (3-2n)/4 , (3+2n)/4 ,  (5-2n)/4 , (5+2n)/4 ; 3/2 ; -1/x² )
 

REGISTERS:    R00 thru R09
FLAGS:            /
 
 

STACK	INPUTS	OUTPUTS
Y	n	Yn(x)
X	x > 0	Jn(x)

     PI     ENTER^
   11.6   XEQ "AEJYNX"  >>>>    J_Pi(11.6)  =  0.238578119
                                        X<>Y   Y_Pi(11.6)  = 0.002890137

-Infinite loop if n = PI and x = 11.5
 

MODIFIED BESSEL FUNCTION 2nd KIND
 

   K_n(x) ~  (pi/(2x))^1/2 e^-x₂F₀( n+1/2 , -n+1/2 ;; -1/(2x) )   for large x-values
 

REGISTERS:    R00 thru R03
FLAGS:            /
 
 

STACK	INPUTS	OUTPUTS
Y	n	/
X	x	Kn(x)

     PI    ENTER^
   10.1  XEQ "AEKNX"  >>>>  K_Pi(10.1)  = 0.00002545492107

-Infinite loop if n = PI and x = 10
 

PARABOLIC CYLINDER FUNCTION
 

    D_n(x)  ~  xⁿ exp(-x²/4) ₂F₀ [ (1-n)/2 , -n/2 ;; -2/x² ]
 

REGISTERS:    R00 thru R04
FLAGS:            /
 
 

STACK	INPUTS	OUTPUTS
Y	n	/
X	x	Dn(x)

    PI     ENTER^
   4.7    XEQ "AEDNX"  >>>>     D_Pi(4.7)  = 0.437982402                     Infinite loop with n = PI & x = 4.6

    PI    CHS
   8.1    XEQ "AEDNX"  >>>>     D_-Pi(8.1)  = 9.594954956 E-11           Infinite loop with n = -PI & x = 8

Note:

-This program may also be used for small x-values when n is a non-negative integer, for instance:

     0     ENTER^
   0.4   XEQ "AEDNX"  >>>>     D₀(0.4)  = 0.960789439
 

KELVIN FUNCTION 1st KIND & 2nd KIND
 

Formulae:     if  x > 0

  ber_n(x) = exp(x/sqrt(2)) (2.PI.x) ^-1/2 [ f_n(x) cos A + g_n(x) sin A ] - (1/PI) [ sin (360°.n) ker_n(x) + cos (360°.n) kei_n(x) ]
  bei_n(x) = exp(x/sqrt(2)) (2.PI.x) ^-1/2 [ f_n(x) sin A - g_n(x) cos A ] + (1/PI) [ cos (360°.n) ker_n(x) - sin (360°.n) kei_n(x) ]

  ker_n(x) = exp(-x/sqrt(2)) (PI/(2x))^1/2 [ f_n(-x) cos B - g_n(-x) sin B ]       with    A = (180°/PI) x/sqrt(2) + ( n/2 - 1/8 ) 180°
  kei_n(x) = exp(-x/sqrt(2)) (PI/(2x))^1/2 [ -f_n(-x) sin B - g_n(-x) cos B ]                B = A + 45°

 and     f_n(x)  ~  1 + Sum_k=1,2,....  [ ( 4n² - 1 ) ( 4n² - 9 ) ...... ( 4n² - ( 2k-1)² ) ] / [ k! (-8x)^k ]   cos (45°k)
           f_n(-x)  ~  1 + Sum_k=1,2,....  [ ( 4n² - 1 ) ( 4n² - 9 ) ...... ( 4n² - ( 2k-1)² ) ] / [ k! (8x)^k ]   cos (45°k)

           g_n(x)  ~   Sum_k=1,2,....  [ ( 4n² - 1 ) ( 4n² - 9 ) ...... ( 4n² - ( 2k-1)² ) ] / [ k! (-8x)^k ]   sin (45°k)
           g_n(-x)  ~  Sum_k=1,2,....  [ ( 4n² - 1 ) ( 4n² - 9 ) ...... ( 4n² - ( 2k-1)² ) ] / [ k! (8x)^k ]   sin (45°k)
 

REGISTERS:    R00 thru R11
FLAGS:           /
 
 

STACK	INPUTS	OUTPUTS
T	/	kein(x)
Z	/	kern(x)
Y	n	bein(x)
X	x > 0	bern(x)

     PI     ENTER^
   10.1   XEQ "AEKLV"       >>>>    ber_Pi(10.1) =  97.94779521
                                             RDN    bei_Pi(10.1) = -56.48175125
                                             RDN    ker_Pi(10.1) =   0.0004270886844
                                             RDN    kei_Pi(10.1) =  -0.00009333763963

-Infinite loop when n = PI & x = 10
 

STRUVE FUNCTION
 

   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² )
 

REGISTERS:    R00 thru R10
FLAGS:            /
 
 

STACK	INPUTS	OUTPUTS
Y	n	/
X	x	Hn(x)

     PI     ENTER^
   16.5   XEQ "AEHNX"  >>>>    H_Pi(16.5)  =  13.35178333

-Infinite loop when n = PI & x = 16.4
 

References:

[1]   Abramowitz and Stegun - "Handbook of Mathematical Functions" -  Dover Publications -  ISBN  0-486-61272-4
[2]   Press , Vetterling , Teukolsky , Flannery - "Numerical Recipes in C++" - Cambridge University Press - ISBN  0-521-75033-4
[3]   http://functions.wolfram.com/
[4]   http://earth-info.nga.mil/GandG/publications/tr8350.2/wgs84fin.pdf