-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    mFp( a1,a2,....,am ; b1,b2,....,bp ; x ) =  SUMk=0,1,2,.....    [(a1)k(a2)k.....(am)k] / [(b1)k(b2)k.....(bp)k] . xk/k!

              where          (ai)k = ai(ai+1)(ai+2) ...... (ai+k-1)   &  (ai)0 = 1    ,    likewise for  (bj)k   ( Pochhammer symbol )
                                   ~
                and           mFp( a1,a2,....,am ; b1,b2,....,bp ; x ) = SUMk=0,1,2,.....    [(a1)k(a2)k.....(am)k] / [Gam(k+b1) Gam(k+b2).....Gam(k+bp)] . xk/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) =   mFp( a1,a2,....,am ; b1,b2,....,bp ; x )    if   Y = + p      hypergeometric function
                  ~
   and     f(x) =   mFp( a1,a2,....,am ; b1,b2,....,bp ; x )    if   Y = - p     regularized hypergeometric function

 >>>   a1,a2,....,am   are to be stored in registers  R01 ............ Rmm
 and    b1,b2,....,b  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+"  >>>>   3F4( 1 , 4 , 7 ; 2 , 3 , 6 , 5 ; Pi )  =  1.631019643

    3   ENTER^
  -4   ENTER^                        ~
   PI  XEQ "HGF+"   >>>>   3F4( 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+"   >>>>   2F2( 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
 

    RF(x,y,z) =  (1/2) §0+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"  >>>   RF(4,3,2) = 0.584082842
 

CARLSON ELLIPTIC INTEGRAL 1ST KIND, CONJUGATE COMPLEX ARGUMENTS
 

  "RFZ"  calculates  RF(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"  >>>   RF(4,3+2.i,3-2.i) = 0.529597761
 

CARLSON ELLIPTIC INTEGRAL 3RD KIND, REAL ARGUMENTS
 

    RJ(x,y,z,p) =  (3/2)  §0+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"  >>>   RJ(4,3,2,1) = 0.360378094

-Alpha register is cleared.
 

CARLSON ELLIPTIC INTEGRAL 3RD KIND, CONJUGATE COMPLEX ARGUMENTS
 

    "RJZ" computes  RJ(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"  >>>   RJ(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
 

     RG(x,y,z) =  (1/4)  §0+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"  >>>   RG(4,3,2) = 1.725503029
 

LEGENDRE ELLIPTIC INTEGRAL 1ST KIND
 

    F ( x | m ) =  §0x ( 1 - m sin2 t ) -1/2 dt  = sin x . RF ( cos2 x ; 1 - m sin2 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 ) =  §0x ( 1 - m sin2 t ) 1/2 dt  = sin x . RF ( cos2 x ; 1 - m.sin2 x ; 1 ) - (m/3) sin3 x . RJ ( cos2 x ; 1 - m.sin2 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 ) =  §0x  ( 1 + n sin2 t ) -1 ( 1 - m sin2 t ) -1/2 dt
                         = sin x . RF ( cos2 x ; 1 - m.sin2 x ; 1 ) - (n/3) sin3 x . RJ ( cos2 x ; 1 - m.sin2 x ; 1 ; 1 + n.sin2 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 =  §0v ( 1 - m sin2 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 sin2 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 RG( a2b2 , a2c2 , b2c2 )
 

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.Pn(x) = (2n-1).x.Pn-1(x) - (n-1).Pn-2(x)  ;  P0(x) = 1  ;  P1(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     P7(4.9) =  1698444.019
                               RDN    P6(4.9) = 188641.3852  is in Y-register.
 

HERMITE POLYNOMIALS
 

     Formulae:      Hn(x) = 2x.Hn-1(x) - 2(n-1).Hn-2(x)  ;  H0(x) = 1  ;  H1(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"  >>>   H7 (3.14) =   73726.24332
                                RDN   H6 (3.14) =   21659.28040
 

CHEBYSHEV POLYNOMIALS - 1ST KIND
 

       Formulae:    Tn(x) = 2x.Tn-1(x) - Tn-2(x)  ;  T0(x) = 1  ;  T1(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"  >>>>   T7(0.314) =  -0.786900700
                                  RDN    T6 (0.314) =    0.338782777
 

CHEBYSHEV POLYNOMIALS - 2ND KIND
 

       Formulae:    Un(x) = 2x.Un-1(x) - Un-2(x)  ;  U0(x) = 1  ;  U1(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"     >>>>    U7 (0.314) =  -0.582815681
                                         RDN    U6 (0.314) =   0.649952293
 

GENERALIZED LAGUERRE POLYNOMIALS
 

        Formulae:      L0(a) (x) = 1  ;   L1(a) (x) = a+1-x   ;    n.Ln(a) (x) = (2.n+a-1-x).Ln-1(a) (x) - (n+a-1).Ln-2(a) (x)
 

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"  >>>>   L7(1.4)(Pi) = 1.688893513
                                 RDN    L6(1.4)(Pi) = 2.271353727
 

ASSOCIATED LEGENDRE FUNCTIONS 1st KIND - INTEGER ORDER & DEGREE
 

       Formulae:      (n-m) Pnm(x) = x (2n-1) Pn-1m(x) - (n+m-1) Pn-2m(x)

                                Pmm(x) = (-1)m (2m-1)!!  ( 1-x2 )m/2   if  | x | < 1                   where (2m-1)!! = (2m-1)(2m-3)(2m-5).......5.3.1
                                Pmm(x) =   (2m-1)!!  ( x2-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"  >>>>     P74(3) =  37920960
                             X<>Y    P64(3) =  2963520

     3   ENTER^
  100  ENTER^
   0.7  XEQ "PMN"  >>>>   P1003(0.7) = -58239.28386
                                X<>Y    P993(0.7) =  13003.91702

-Actually,  Pnm(x)  is not always a polynomial !
 

JACOBI POLYNOMIALS
 

       Formulae:      P0(a;b) (x) = 1  ;   P1(a;b) (x) = (a-b)/2 + x.(a+b+2)/2

                             2n(n+a+b)(2n+a+b-2).Pn(a;b) (x) = [ (2n+a+b-1).(a2-b2) + x.(2n+a+b-2)(2n+a+b-1)(2n+a+b) ] Pn-1(a;b) (x)
                                                                                    - 2(n+a-1)(n+b-1)(2n+a+b) Pn-2(a;b) (x)
 

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"  >>>>   P7(1.4;1.7)(1/Pi) = -0.322234421
                                           X<>Y  P6(1.4;1.7)(1/Pi) =   0.538220923
 

ULTRASPHERICAL POLYNOMIALS
 

           Formulae:      C0(a) (x) = 1  ;   C1(a) (x) = 2.a.x   ;    (n+1).Cn+1(a) (x) = 2.(n+a).x.Cn(a) (x) - (n+2a-1).Cn-1(a) (x)      if  a # 0

                                  Cn(0) (x) = (2/n).Tn(x)
 

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 )

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

Examples:

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

   0   ENTER^
   7   PI   1/X   XEQ "USP"  >>>>   C7(0)(1/Pi)  =  -0.219109374
 

-MISCELLA
 

LINEAR PROGRAMMING: SIMPLEX METHOD
 

   SIMPLEX solves linear programming problems:

The purpose is to find  m  non-negative real numbers:   x1 , ..... , xm
satisfying:
                          bi ³ ai;1x1 + ....... + ai;m xm            i = 1 , .... , n                                ( n inequations )                   all the bi  , bi'  and  bi"
                          bi' = ai';1x1 + ....... +ai';m xm            i' = 1 , ..... , n'                             ( n' equations )                    must be non-negative*
                          bi"£ ai";1x1+ ....... +ai";m xm         i" = 1 , ..... , n"                             (n" inequations )

and such that      F  =  c1x1 + ....... +  cmxm       is maximized.    ( to find the minimum of F, seek the maximum of  -F )
 

 *  if some bi 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 cj  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 §x1xn 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 = x1       R03 = x2      .......    R2n-1 = xn
       R02 = y1       R04 = y2      .......     R2n   = yn

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

Example:     Calculate  §18  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"    >>>>    §18  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:   M1 ( x1 , y1 )  ,   M2 ( x2 , y2 )  ,  ..........  ,  Mn ( xn , yn )
  store the 2n numbers   x1 , y1 , x2 , y2 , ........ , xn , yn  into contiguous registers,  starting at any register  Rbb:

       Rbb = x1         Rb+2 = x2      .......       Ree-1 = xn
       Rb+1 = y1       Rb+3 = y2      .......       Ree   = yn                  ( 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 §x1x2 L(x) dx , given a set of n data points which must be stored as follows:

       R00 = n = number of points

       R01 = x1       R03 = x2      .......    R2n-1 = xn
       R02 = y1       R04 = y2      .......     R2n   = yn

REGISTERS:    R00 thru R2n
FLAGS:            /
 
 
      STACK        INPUT      OUTPUT
           X              /    §x1xn  L(x).dx

Example:     Calculate  §18  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"     >>>>    §18  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(xi,yj)   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/103

                 x  \  y |         Rnn+1 = y1               R2n+2 = y2       ..............................      Rnm+m = ym
         ---------------------------------------------------------------------------------------------------------------------
          R01  = x1 |     Rnn+2 = f(x1,y1)      R2n+3 = f(x1,y2)    ..............................       Rnm+m+1 = f(x1,ym)
          R02  = x2 |     Rnn+3 = f(x2,y1)      R2n+4 = f(x2,y2)    ..............................       Rnm+m+2 = f(x2,ym)
             ...........    |      .............................................................................................................................................
                            |
           Rnn = xn  |     R2n+1  = f(xn,y1)     R3n+2 = f(xn,y2)   ..............................        Rnm+m+n  = f(xn,ym)
 

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 Snm are defined by the recurrence relation: ( Sn0 = 0 ) ;  Snm = Sn-1m-1 - ( n-1 ) Sn-1m      ( 1 <= m <= n )
-Stirling numbers of the second kind snm ----------------------------------  ( sn0 = 0 )  ;   snm =  sn-1m-1  +  m sn-1m            --------------
 

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  S127   and    s127

   a)  CF 02
        12  ENTER^
          7  XEQ "SNK"   yields       S127 =  -2637558      ( and R01 = -39916800 = S121 ; R02 = 120543840 = S122
                                                                                          .........................................      R07 = -2637558 = S127  )

   b)  SF 02
        12  ENTER^
         7     R/S          produces       s127 = 627396       ( and R01 = 1 =  s121 ; R02 = 2047 =  s122 ; ........ ;  R07 = 627396 = s127 )
 
 

CLEBSH-GORDAN COEFFICIENTS & WIGNER 3J-SYMBOL
 

     Clebsch-Gordan Coefficients  = <  j1 j2  m1 m2 |   j1 j2 j m  >
                                                    = [ ( 2j +1) D( j1 j2 j  ) ]1/2  [ ( j1 + m1 )! ( j1 - m1 )! ( j2 + m2 )! ( j2 - m2 )! ( j + m )! ( j - m )! ]
                                   x Sumk  (-1)k / [ k! ( j1 + j2 - j - k )! ( j1 - m1 - k )! ( j2 + m2 - k )! ( j - j2 + m1 + k )! ( j - j1 - m2 + k )! ]         if   m1 + m2 = m

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

     <  j1 j2  m1 m2 |   j1 j2 j m  > = 0            if   m1 + m2 # m

     Wigner 3j-symbol  = W3J (  j1 j2 j  m1 m2 -m ) = (-1)^( j1 - j2 + m )  ( 2j +1) -1/2  <  j1 j2  m1 m2 |   j1 j2 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 = j1       R02 = j2      R03 = j
       R04 = m1     R05 = m2     R06 = m

REGISTERS:    R00 thru R14
FLAGS:            F10
 
 
      STACK        INPUTS      OUTPUTS
           Y             /          W3j 
           X             /          ClG

   where    W3j = W3j (  j1j2 j  m1 m2 -m )      and  ClG = <  j1 j2  m1 m2 |   j1 j2 j m  >

Examples:

-If   j1 = 12  ,   j2 = 24  ,  j  = 31
      m1 = 1   ,  m2 = 16  , m = 17

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

-With   j1 = 1  ,   j2 = 3/2  ,  j  = 5/2
           m1 = 0 ,  m2 = 3/2  , m = 3/2

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

WIGNER 6J-SYMBOL
 

     W6J ( j1 j2 j3 j4 j5 j6 )  =  [ D ( j1 j2 j3 ) D ( j1 j5 j6 ) D ( j2 j3 j4 ) D ( j3 j4 j5 ) ]1/2
                                       x Sumk=k' to k"  ( -1 )k ( k + 1 )! / [ ( k - j1 - j2 - j3 )! ( k - j1 - j5 - j6 )! ( k - j2 - j3 - j4 )! ( k - j3 - j4 - j5 )!
                                                                                              ( j1 + j2 + j4 + j5 - k )! ( j1 + j3 + j4 + j6 - k )! ( j2 + j3 + j5 + j6 - k )! ]

       with   k' = max { j1 + j2 + j3 , j1 + j5 + j6 , j2 + j3 + j4 , j3 + j4 + j5 }
       and    k" = min { j1 + j2 + j4 + j5 , j1 + j3 + j4 + j6 , j2 + j3 + j5 + j6 }
 

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

          R01 = j1       R02 = j2      R03 = j3
          R04 = j4       R05 = j5      R06 = j6

REGISTERS:    R00 thru R14
FLAGS:            F10
 
 
      STACK        INPUTS      OUTPUTS
           X             /          W6j

    with  W6j = W6j ( j1 j2 j3 j4 j5 j6 )

Example:

-If   j1 = 10  ,   j2 = 16  ,  j3  = 21   which are to be stored into      R01   R02   R03
       j4 = 24  ,   j5 = 12  ,  j6  = 14                                                 R04   R05   R06

 XEQ "W6J"  >>>>   W6j  =  0.006251585
 

WIGNER 9J-SYMBOL
 

     W9J ( j1 j2 j3 j4 j5 j6 j7 j8 j9 )  =  Sumk=k' to k"  (-1)2k (2k+1)  W6J ( j1 j9 k j8 j4 j7 )  W6J ( j8 j4 k j6 j2 j5 )  W6J ( j6 j2 k j1 j9 j3 )
 

    where    k' = max { | j1 - j9 | , | j4 - j8 | , | j2 - j6 | }    and    k" =  min {  j1 + j9  ,  j4 + j8  ,  j2 + j6 }
 

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

         R01 = j1       R02 = j2      R03 = j3
         R04 = j4       R05 = j5      R06 = j6              registers R01 to R09 are saved in registers R15 to R23
         R07 = j7       R08 = j8      R09 = j9              and they are restored at the end

REGISTERS:    R00 thru R27
FLAGS:            F10
 
 
      STACK        INPUTS      OUTPUTS
           X             /          W9j

  where  W9j = W9J ( j1 j2 j3 j4 j5 j6 j7 j8 j9 )

Examples:

         j1 = 3  ,   j2 = 7  ,  j3  = 5                                                   R01   R02   R03
-If     j4 = 6  ,   j5 = 8  ,  j6  = 9     which are to be stored into      R04   R05   R06
         j7 = 4  ,   j8 = 5  ,  j9  = 7                                                   R07   R08   R09

 XEQ "W9J"  >>>>   W9j  =  0.001032402067

               j1 = 1/2  ,   j2 =  1    ,  j3  = 3/2
-With      j4 =  1    ,   j5 = 1/2  ,  j6  = 1/2
               j7 = 3/2  ,   j8 = 1/2  ,  j9  =  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:   §ab §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 = §12 §xx^2 (1+x4 y4)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:      §ab §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 = §12 §xx^2 §x+yxy  xyz  ( x2 + y2 + z2 ) -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    I1 =  0.765014888
         n = 2      2  STO 03         R/S            "        I2 =  0.770640690
         n = 4      4  STO 03         R/S            "        I4 =  0.770731245
         n = 8      8  STO 03         R/S            "        I8 =  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 = ( L2 - L1 )/2  ;   F = ( b1 + b2 )/2  ;  G = ( b2 - b1 )/2
           S = sin2 G  cos2 L  +  cos2 F  sin2 L  ;   µ = Arcsin S1/2  ;  R = ( S.(1-S) )1/2

 we define  A = ( sin2 G  cos2 F )/S + ( sin2 F  cos2 G )/( 1-S )
         and   B = ( sin2 G  cos2 F )/S - ( sin2 F  cos2 G )/( 1-S )

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

  where   eps =    (f/2) ( -A - 3B R/µ )
                       + (f2/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 Az1 in the direction from P1 towards P2 is obtained by the following formulae:    Az1 = Az1sph + dAz

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

       and    Tan Az1sph = cos b2 sin 2L / ( cos b1 sin b2 - sin b1 cos b2 cos 2L )     >>>  the R-P function returns the angle in the proper quadrant

  Az1sph 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 P2 to P1 ) is obtained by replacing L by -L and exchanging b1 & b2 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    Az1 =  51°47'36"76 = forward azimuth
                                          RDN    Az2 = -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  =  ( p2+ q2 )1/2

 where   p =  { - [ GM + ( omega )2 a2 E ( q' / 2q0 ) ( sin2 F - 1/3 ) ] / ( u2 + E2 ) + ( omega )2 u cos2 F } / w
  and     q =   [ a2 ( q / q0 ) ( u2 + E2 ) -1/2 - ( u2 + E2 )1/2 ] ( omega )2 ( sin F cos F ) / w

            u  =  { [ ( x2 + y2 + z2 - E2 ) / 2 ] [ 1 + { 1 + 4 E2 z2 / ( x2 + y2 + z2 - E2 )2 }1/2 ] }1/2
  with    E =   a ( 2 f - f 2 )1/2  ,   sin F = z / u   ,   cos F = ( x2 + y2 )1/2 / ( u2 + E2 )1/2
            w =  [ ( u2 + E2 sin2 F )1/2 / ( u2 + E2 )1/2 ]1/2

            q  = (1/2) [ ( 1 + 3 u2 / E2 ) Atan ( E / u )  - 3 E / u ]                                                            ( x2 + y2 )1/2  = ( N + h ) cos b
  and    q0 = (1/2) [ ( 1 + 3 b' 2 / E2 ) Atan ( E / b' ) - 3 E / b' ]         ( b' = polar radius )                          z            = [ ( b'/a )2 N + h ] sin b
            q' =   3 ( 1 + u2 / E2 ) [ 1 - ( u / E ) Atan ( E / u ) ] - 1                                                                 N          = a ( 1 - e2 sin2 b ) -1/2
 

-If, however, we use ATAN to calculate  q , q0 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 t3  Sumk=0,1,2,.....  ( k + 1 ) ( - t2 )k / [ ( 2 k + 3 ) ( 2 k + 5 ) ]
    3 ( 1 + 1 / t2 ) [ 1 - ( Atan t ) / t ] - 1  =  6 t2  Sumk=0,1,2,.....  ( - t2 )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/s2 ) 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/s2
                                         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/s2
                                    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) ] 0F1( ; 2/3 ; x3/9 ) - x [ 3 -1/3 / Gam(1/3) ] 0F1( ; 4/3 ; x3/9 )
    Bi(x) = [ 3 -1/6 / Gam(2/3) ] 0F1( ; 2/3 ; x3/9 ) + x [ 3 1/6 / Gam(1/3) ] 0F1( ; 4/3 ; x3/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/pi1/2)  §0x  e-t^2 dt

      Ascending series                 erf x = (2x/Pi1/2) exp(-x21F1( 1 , 3/2 , x2 )

    Continued Fraction:      2.ex^2  §xinfinity 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) = §0x  sin(pi.t2/2).dt      and      C(x) = §0x  cos(pi.t2/2).dt

Formulae:       S(x) = ( PI x3/6 ) 1F2( 3/4 ; 3/2 , 7/4 ; -PI2 x4/16 )

                        C(x) = x 1F2( 1/4 ; 1/2 , 5/4 ; -PI2 x4/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 z2 =  (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)  = §0x sin(t)/t dt = x 1F2( 1/2 ; 3/2 , 3/2 ; -x2/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) = §0x sinh(t)/t dt = x 1F2( 1/2 ; 3/2 , 3/2 ; x2/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) + §0x (cos(t)-1)/t dt = C + Ln(x) - (x2/4) 2F3( 1 , 1 ; 2 , 2 , 3/2 ; -x2/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) + §0x (cosh(t)-1)/t dt = C + Ln(x) + (x2/4) 2F3( 1 , 1 ; 2 , 2 , 3/2 ; x2/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-12/(3+i.x-22/(5+i.x-32/(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) = §-infinityx  et/t dt  = C + Ln(x) + x 2F2( 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 En(x)
 

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

     En(x) = (-x)n-1 ( -Ln x - C + Sumk=1,...,n-1 1/k ) / (n-1)!  - Sumk#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

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

-We also have:   E0(x) = (1/x).exp(-x)
  and  En(x) = xn-1 Gam(1-n) - [1/(1-n)] 1F1( 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"  >>>>  E0(1.4) = 0.176140689

       2    ENTER^
     1.6  XEQ "ENX"  >>>>  E2(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/2x + 1/3x + ........ + 1/nx + ......     ( for x > 1)
  but the following program uses the formula:          Zeta(x) = ( 1 - 1/2x + 1/3x - 1/4x + ...... ) / ( 1 - 21-x )   where roundoff errors are smaller.

-If  x < 1, we have:   Zeta(x) = GAM(1-x)  2x  PIx-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

                             ~
    In(x) = ( x/2 )n 0F1( ; n+1 ; x2/4 )
 

REGISTERS:    R00 thru R03
FLAGS:            F01
 
 
      STACK        INPUTS     OUTPUTS
           Y             n            /
           X             x          In(x)

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

-The results are accurate for large arguments too.
 

BESSEL FUNCTION 1st KIND

                            ~
   Jn(x) = ( x/2 )n 0F1( ; n+1 ; - x2/4 )
 

REGISTERS:    R00 thru R03
FLAGS:            F01
 
 
      STACK        INPUTS     OUTPUTS
           Y             n            /
           X             x          Jn(x)

   0.7   ENTER^
   1.9   XEQ "JNX"   >>>>   J0.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:   J0(x) + 2 J2(x) + 2 J4(x) + 2 J6(x) + ...... = 1  and   Jn-1(x) + Jn+1(x) = (2n/x) Jn(x)
 

REGISTERS:    R00 thru R05
FLAGS:            /
 
 
      STACK        INPUTS     OUTPUTS
           Y             n            /
           X             x          Jn(x)

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

BESSEL FUNCTION 2nd KIND
 

  Yn(x) = ( Jn(x) cos(n(pi)) - J-n(x) ) / sin(n(pi))    if n is not an integer and otherwise ( n non-negative ):

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

REGISTERS:    R00 thru R05
FLAGS:            F01
 
 
      STACK        INPUTS     OUTPUTS
           Y             n            /
           X             x          Yn(x)

Examples

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

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

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

MODIFIED BESSEL FUNCTION 2nd KIND
 

      Kn(x) = (pi/2) ( I-n(x) - In(x) ) / sin(n(pi))      if n is not an integer and otherwise ( n non-negative ):

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

REGISTERS:    R00 thru R05
FLAGS:            F01
 
 
      STACK        INPUTS     OUTPUTS
           Y             n            /
           X             x          Kn(x)

Examples

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

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

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

PARABOLIC CYLINDER FUNCTION
 

   Dn(x) = 2n/2 Pi1/2 exp(-x2/4) [ 1/Gam((1-n)/2) 1F1( -n/2 , 1/2 , x2/2 ) - 21/2 ( x / Gam(-n/2) ) 1F1[ (1-n)/2 , 3/2 , x2/2 ]
 

REGISTERS:    R00 thru R05
FLAGS:            /
 
 
      STACK        INPUTS      OUTPUTS
           Y             n             /
           X             x         Dn(x)

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

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

KELVIN FUNCTION 1st KIND
 

Formulae:

   bern(x) = [ cos(135°n) / Gam(n+1) ] (x/2)n S  -  [ sin(135°n) / Gam(n+2) ] (x/2)n+2 T

   bein(x) = [ cos(135°n) / Gam(n+2) ] (x/2)n+2 S  +  [ sin(135°n) / Gam(n+1) ] (x/2)n T
 

 where      S = 0F3( ; (n+1)/2 , 1+n/2 ,1/2 ; - x4/256 )       T = 0F3( ; 1+n/2 , (n+3)/2 , 3/2 ; - x4/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"  >>>>    bersqrt(2)(PI) = -0.674095953
                            X<>Y    beisqrt(2)(PI) = -1.597357212

-For large arguments, "AEKLV" is preferable.
 

KELVIN FUNCTION 2nd KIND
 

-If n is not an integer, we have:

     kern(x) = - (PI/2) [ bein(x) - (ber-n(x)/sin(n.180°)) + (bern(x)/tan(n.180°)) ]

    kein(x)  =   (PI/2) [ bern(x) + (bei-n(x)/sin(n.180°)) - (bein(x)/tan(n.180°)) ]

-Otherwise,

     kern(x) = (1/2) (x/2)-n SUMk=0,1,....,n-1 cos ((3n/4+k/2)180°) (n-k-1)!/(k!) (x2/4)k - ln(x/2) bern(x) + (PI/4) bein(x)
                 + (1/2) (x/2)n SUMk=0,1,.....   cos ((3n/4+k/2)180°)  ( psi(k+1) + psi(n+k+1) ) (x2/4)k / (k!(n+k)!)

     kein(x) = - (1/2) (x/2)-n SUMk=0,1,....,n-1 sin ((3n/4+k/2)180°) (n-k-1)!/(k!) (x2/4)k - ln(x/2) bein(x) - (PI/4) bern(x)
                   + (1/2) (x/2)n SUMk=0,1,.....   sin ((3n/4+k/2)180°)  ( psi(k+1) + psi(n+k+1) ) (x2/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"  >>>>    kersqrt(2)(PI) =  0.025901894
                               X<>Y    keisqrt(2)(PI) =  0.089242866

   3    ENTER^
  PI    XEQ "KLV2"    >>>>    ker3(PI) =  -0.045754846
                                  X<>Y    kei3(PI) =  -0.196928536

-For large arguments, "AEKLV" is preferable.
 

STRUVE FUNCTION

                                 ~
    Hn(x) =  (x/2)n+1 1F2( 1 ; 3/2 , n+3/2 ; - x2/4 )
 

REGISTERS:    R00 thru R04
FLAGS:            F01
 
 
      STACK        INPUTS     OUTPUTS
           Y             n            /
           X             x          Hn(x)

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

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

MODIFIED STRUVE FUNCTION

                                 ~
   Ln(x) =    (x/2)n+11F2( 1 ; 3/2 , n+3/2 ; x2/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   Bx(a,b) = §0x  ta-1 (1-t)b-1 dt       ( a , b > 0 )
-We have also  Bx(a,b) = (xa/a) (1-x)b2F1(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"  >>>>  B0.7(Pi,e) =  0.029623046
 

INCOMPLETE GAMMA FUNCTION
 

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

     igam(a,x) =  (xa/a).exp(-x) 1F1(1,a+1,x)    where    1F1 = 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     Jn(x) = (1/PI) §0pi  cos ( n.t - x.sin t ) dt
-Weber's function is defined by    En(x) = (1/PI) §0pi  sin ( n.t - x.sin t ) dt

-The following programs use the formulae:
                                               ~                                                                         ~
   Jn(x) = + (x/2) sin ( 90°n )  1F2( 1 ; (3-n)/2 , (3+n)/2 ; -x2/4 )  + cos ( 90°n ) 1F2( 1 ; (2-n)/2 , (2+n)/2 ; -x2/4 )
                                                ~                                                                       ~
   En(x) = - (x/2) cos ( 90°n )  1F2( 1 ; (3-n)/2 , (3+n)/2 ; -x2/4 ) + sin ( 90°n ) 1F2( 1 ; (2-n)/2 , (2+n)/2 ; -x2/4 )
 

REGISTERS:    R00 thru R06
FLAGS:            /
 
 
      STACK        INPUTS      OUTPUTS
           Y             n         Jn(x)
           X             x         En(x)

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

DEBYE FUNCTION
 

   dby(x;n) = §x+infinity  tn/(et-1).dt

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

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

                                            ~
    Pnm(x) = |(x+1)/(x-1)|m/22F1(-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"  >>>>   P1.30.4(0.7) =  0.274932821

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

ASSOCIATED LEGENDRE FUNCTIONS 2nd KIND
 

Formulae:

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

   Qnm(x) = ei.mpi 2-n-1 pi1/2 x-n-m-1 (x2-1)m/2 2F1(1+n/2+m/2 ; 1/2+n/2+m/2 ; n+3/2 ; 1/x2 )  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"   >>>>     Q1.30.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,     Q1.20.7(1.9) = -0.081513570 + 0.112193805 i
 

F-DISTRIBUTION FUNCTION
 

  FDF   calculates    Q( x | n1 , n2 ) =  (n1/n2)(n1)/2  Gamma[(n1+n2)/2] / [ Gamma(n1/2) Gamma(n2/2) ] §xinfinity  t -1+(n1)/2 ( 1 + t.n1/n2 ) -(n1+n2)/2 dt

  via the incomplete beta function:   Q( x | n1 , n2 ) = Iy( n2/2 , n1/2 )   with  y = n2/(n2+n1.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( X2 | m ) = 2 -m/2 [ 1/Gamma(m/2) ] §0X^2   t -1+m/2 e -t/2 dt

-It uses the series exansion:    P( X2 | m ) = (X2/2)m/2 e -(X^2)/2 [ 1/Gamma(m/2) ]  [ 1 + Sumk=1,2,...  X2k/[ (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 + x2/n ) -(n+1)/2  dx

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

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

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:         FL(n,r) = CL(n) r L+1 Sum k>L AkL (n) r k-L-1

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

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

    | Gam( 1+L+i y ) |2 = [ L2 + y2 ] [ (L-1)2 + y2 ] .................. [ 1 + y2 ]  (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"  >>>>   F2( 0.7 , 1.8 ) = 0.141767746
 

TORONTO FUNCTION
 

Toronto function is defined by

   T(m,n,r) = exp(-r2) [ Gam((m+1)/2) / n! ] r2n-m+11F1( (m+1)/2 , n+1 , r2 )
 

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 Lmn by solving the transcendental equation   U1(Lmn) + U2(Lmn) = 0   where    U1 & U2  are 2 continued fractions:

                                                 -brm                      -br-2m
   U1(Lmn) = grm - Lmn +    ---------------       ---------------   ..............
                                           gr-2m - Lmn +         gr-4m - Lmn +

                                                                                                                  r = n - m
                                     -br+2m                    -br+4m
   U2(Lmn) =           ---------------       ----------------   ..............
                               gr+2m - Lmn +         gr+4m - Lmn +
 

   with      brm  =  [ r.(r-1).(2m+r).(2m+r-1).c4 ] / [ (2m+2r-1)2.(2m+2r+1).(2m+2r-3) ]
   and      grm  =  (m+r).(m+r+1) + (c2/2).[ 1 - (4m2-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 ,  c2 = -1    ( note that  c2  may be positive or negative )

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

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

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

Notes:

-Usually, m & n are non-negative integers and n = m , m+ 1 , m + 2 , ........
  c2 = 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 , c2 , 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
 

  Smn(x) = ( 1 - x2 )m/2  f(x)        where   f(x) = Sumk=0,1,.... ak xk

     with  a0 = Pnm(0)  =  2m sqrt(PI) / [ Gam((1-m-n)/2) Gam((2-m+n)/2 ]                           ( Flammer's scheme )
     and   a1 = P'nm(0) = ( m + n ) 2m sqrt(PI) / [ Gam((2-m-n)/2) Gam((1-m+n)/2 ]

-In both cases,   (k+1)(k+2) ak+2 - [ k ( k + 2m + 1 ) - Lmn + m ( m + 1 ) ] ak  - c2  ak-2 = 0
 

REGISTERS:    R00 thru R10
FLAGS:           F10
 
 
      STACK        INPUT      OUTPUT
           X             x         Smn(x)

  Provided that m , n , c2 , 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   c2 = 1.7    "LMNC2"  gives   2.246866650       -Store these 4 numbers into  R01 thru R04  and

    0.7   XEQ "SWF"  >>>>     Smn(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  2F1(a,b,c,x) = 1 +  (a)1(b)1/(c)1. x1/1! + ............. +  (a)n(b)n/(c)n . xn/n! + ..........  where (a)n = a(a+1)(a+2) ...... (a+n-1)

-The relation        2F1(a,b,c,x) =  ( 1 - x ) -a  2F1(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  2F1(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 x3/22F0( 1/6 , 5/6 ; -3/(4x3/2) )
     Bi(x) ~ { 1/sqrt[ pi sqrt(x) ]  }  exp ( 2/3 x3/22F0( 1/6 , 5/6 ; 3/(4x3/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) x3/2    ;     S = Sumk=0,1,2,.... (-1)k c2k z -2k    ,    T = Sumk=0,1,2,.... (-1)k c2k+1 z -2k-1

       and       c0 = 1 ,  ck = ck-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:

     Jn(x) ~ (2/(pi*x))1/2 ( P.cos(x-(2n+1)pi/4) + Q.sin(x-(2n+1)pi/4) )
    Yn(x) ~ (2/(pi*x))1/2 ( P.sin(x-(2n+1)pi/4) - Q.cos(x-(2n+1)pi/4) )

  where   P =  4F1( (1-2n)/4 , (1+2n)/4 ,  (3-2n)/4 , (3+2n)/4 ; 1/2 ; -1/x2 )
     and   Q =  [(1-4n2)/(8x)] 4F1( (3-2n)/4 , (3+2n)/4 ,  (5-2n)/4 , (5+2n)/4 ; 3/2 ; -1/x2 )
 

REGISTERS:    R00 thru R09
FLAGS:            /
 
 
      STACK        INPUTS     OUTPUTS
           Y             n          Yn(x)
           X          x > 0          Jn(x)

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

-Infinite loop if n = PI and x = 11.5
 

MODIFIED BESSEL FUNCTION 2nd KIND
 

   Kn(x) ~  (pi/(2x))1/2 e-x2F0( 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"  >>>>  KPi(10.1)  = 0.00002545492107

-Infinite loop if n = PI and x = 10
 

PARABOLIC CYLINDER FUNCTION
 

    Dn(x)  ~  xn exp(-x2/4) 2F0 [ (1-n)/2 , -n/2 ;; -2/x2 ]
 

REGISTERS:    R00 thru R04
FLAGS:            /
 
 
      STACK        INPUTS      OUTPUTS
           Y             n             /
           X             x         Dn(x)

    PI     ENTER^
   4.7    XEQ "AEDNX"  >>>>     DPi(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"  >>>>     D0(0.4)  = 0.960789439
 

KELVIN FUNCTION 1st KIND & 2nd KIND
 

Formulae:     if  x > 0

  bern(x) = exp(x/sqrt(2)) (2.PI.x) -1/2 [ fn(x) cos A + gn(x) sin A ] - (1/PI) [ sin (360°.n) kern(x) + cos (360°.n) kein(x) ]
  bein(x) = exp(x/sqrt(2)) (2.PI.x) -1/2 [ fn(x) sin A - gn(x) cos A ] + (1/PI) [ cos (360°.n) kern(x) - sin (360°.n) kein(x) ]

  kern(x) = exp(-x/sqrt(2)) (PI/(2x))1/2 [ fn(-x) cos B - gn(-x) sin B ]       with    A = (180°/PI) x/sqrt(2) + ( n/2 - 1/8 ) 180°
  kein(x) = exp(-x/sqrt(2)) (PI/(2x))1/2 [ -fn(-x) sin B - gn(-x) cos B ]                B = A + 45°

 and     fn(x)  ~  1 + Sumk=1,2,....  [ ( 4n2 - 1 ) ( 4n2 - 9 ) ...... ( 4n2 - ( 2k-1)2 ) ] / [ k! (-8x)k ]   cos (45°k)
           fn(-x)  ~  1 + Sumk=1,2,....  [ ( 4n2 - 1 ) ( 4n2 - 9 ) ...... ( 4n2 - ( 2k-1)2 ) ] / [ k! (8x)k ]   cos (45°k)

           gn(x)  ~   Sumk=1,2,....  [ ( 4n2 - 1 ) ( 4n2 - 9 ) ...... ( 4n2 - ( 2k-1)2 ) ] / [ k! (-8x)k ]   sin (45°k)
           gn(-x)  ~  Sumk=1,2,....  [ ( 4n2 - 1 ) ( 4n2 - 9 ) ...... ( 4n2 - ( 2k-1)2 ) ] / [ 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"       >>>>    berPi(10.1) =  97.94779521
                                             RDN    beiPi(10.1) = -56.48175125
                                             RDN    kerPi(10.1) =   0.0004270886844
                                             RDN    keiPi(10.1) =  -0.00009333763963

-Infinite loop when n = PI & x = 10
 

STRUVE FUNCTION
 

   Hn(x)  ~  Yn(x)  + [ 21-n xn-1 / sqrt(Pi) / Gam(n+1/2) ]  3F0( 1/2 , -n + 1/2 , 1 ;; -4/x2 )
 

REGISTERS:    R00 thru R10
FLAGS:            /
 
 
      STACK        INPUTS     OUTPUTS
           Y             n            /
           X             x          Hn(x)

     PI     ENTER^
   16.5   XEQ "AEHNX"  >>>>    HPi(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