Anionic Special Functions(III) for the HP-41

Overview
 

 1°)  Polygamma Functions
 2°)  Incomplete Gamma Function
 3°)  Incomplete Beta Function
 4°)  Harmonic Numbers
 5°)  Bessel-Clifford Functions
 6°)  Whittaker-M Functions
 7°)  Whittaker-W Functions
 8°)  Generalized Error-Functions
 9°)  Lambert W-Function
10°) Continued Fractions ( with an application to the error-function )
 
 

1°)  Polygamma Functions
 

Formulae:     "PSINA" employs the asymptotic expansions:

  •   If m > 0  ,  y^(m) (a)  ~  (-1)^m-1 [ (m-1)! / a^m + m! / (2.a^m+1) + SUM_k=1,2,.... B_2k (2k+m-1)! / (2k)! / a^2k+m       where  B_2k are Bernoulli numbers

  •   If m = 0  ,  y (a)  ~  Ln a - 1/(2a) - SUM_k=1,2,.... B_2k / (2k) / a^2k              ( digamma function )
 

 and the recurrence relation:   y^(m) (a+p) = y^(m) (a) + (-1)^m m! [ 1/a^m+1 + ....... + 1/(a+p-1)^m+1 ]           where  p  is a positive integer
 
 

Data Registers:           •  R00 = n > 1         ( Registers R00 thru Rnn are to be initialized before executing "PSINA" )

                                      •  R01   ......  •  Rnn = the n components of the anion

                                         Rn+1 ..........  R4n:  temp

  >>>  When the program stops:     R01   ......  Rnn = the n components of the result
 

Flags:  /
Subroutines:  "A*A1"  "1/A"  "A^X"   "LNA"    ( cf "Anions for the HP41" )
                         "ST*A"  "ST/A"                           ( cf "Anionic Special Functions(I) for the HP41" )
 
 

 
   ( 401 bytes / SIZE 4n+1 )
 
 

   where  m  is a non-negative integer  <  61  and  1.nnn   is the control number of  the result

Example1:    Calculate  y⁽³⁾ ( 1 + 0.9 i + 0.8 j + 0.7 k )                        ( quaternion ->  4  STO 00 )

     1   STO 01
   0.9  STO 02
   0.8  STO 03
   0.7  STO 04

    3    XEQ "PSINA"   >>>>   1.004                                           ---Execution time = 61s---

   R01 = -0.673649101
   R02 =  0.147195368
   R03 =  0.130840327
   R04 =  0.114485286

  whence    y⁽³⁾ ( 1 + 0.9 i + 0.8 j + 0.7 k ) =  -0.673649101 + 0.147195368 i + 0.130840327 j + 0.114485286 k
 

Example2:    Calculate  y ( 1 + 0.9 i + 0.8 j + 0.7 k )       digamma function    i-e   m = 0              ( quaternion ->  4  STO 00 )

     1   STO 01
   0.9  STO 02
   0.8  STO 03
   0.7  STO 04

    0    XEQ "PSINA"   >>>>   1.004                                           ---Execution time = 47s---

   R01 = 0.377339972
   R02 = 0.783351925
   R03 = 0.696312822
   R04 = 0.609273719

  whence    y ( 1 + 0.9 i + 0.8 j + 0.7 k ) =  0.377339972 + 0.783351925 i + 0.696312822 j + 0.609273719 k
 

Note:

-"PSINA" computes  (m+9)!  - line 66 - so m cannot exceed 60.
 

2°)  Incomplete Gamma Function
 

Formula:      g(m,a) = ( a^m / m ) exp(-a)  M(1,m+1;a)   where  M = Kummer's function
 
 

Data Registers:           •  R00 = n > 1                ( Registers R00 thru Rnn are to be initialized before executing "IGFA" )

                                      •  R01   ......  •  Rnn = the n components of the anion

                                         Rn+1 ..........  R3n:  temp     R3n+1 = m

  >>>  When the program stops:     R01   ......  Rnn = the n components of the result
 

Flags:  /
Subroutines:  "A*A1"  "A^X"   "E^A"              ( cf "Anions for the HP41" )
                         "ST*A"  "ST/A"  "KUMA"          ( cf "Anionic Special Functions(I) for the HP41" )
 
 
 

 
      ( 103 bytes / SIZE 3n+2 )
 
 

   where   1.nnn   is the control number of  the result

Example:    Compute   g( 1.6 , 1 + 2 i + 3 j + 4 k )               ( quaternion ->  4  STO 00 )

     1   STO 01
     2   STO 02
     3   STO 03
     4   STO 04

    1.6   XEQ "IGFA"  >>>>    1.004                                           ---Execution time = 38s---

   R01 =  0.955358734
   R02 = -0.390239936
   R03 = -0.585359904
   R04 = -0.780479873

-So,     g( 1.6 , 1 + 2 i + 3 j + 4 k ) = 0.955358734 - 0.390239936 i - 0.585359904 j - 0.780479873 k
 

3°)  Incomplete Beta Function
 

Formula:    B_a(p,q) = ( a^p / p )  F(p,1-q;p+1;a)   where  F = the hypergeometric function
 

Data Registers:           •  R00 = n > 1                ( Registers R00 thru Rnn are to be initialized before executing "IBFA" )

                                      •  R01   ......  •  Rnn = the n components of the anion

                                         Rn+1 ..........  R3n:  temp

  >>>  When the program stops:     R01   ......  Rnn = the n components of the result
 

Flag:  F04
Subroutines:  "A*A1"  "A^X"                          ( cf "Anions for the HP41" )
                         "ST/A"  "HGFA"          ( cf "Anionic Special Functions(I) for the HP41" )
 
 

 
( 61 bytes / SIZE 3n+1 )
 
 

   where   1.nnn   is the control number of  the result

Example:    Compute  B_{0.1+0.2i+0.3j+0.4k} ( 0.7 , 1.8 )               ( quaternion ->  4  STO 00 )

   0.1  STO 01
   0.2  STO 02
   0.3  STO 03
   0.4  STO 04

   0.7  ENTER^
   1.8  XEQ "IBFA"  >>>>    1.004                                           ---Execution time = 32s---

   R01 = 0.653131938
   R02 = 0.244542459
   R03 = 0.366813688
   R04 = 0.489084918

-Whence,    B_{0.1+0.2i+0.3j+0.4k} ( 0.7 , 1.8 ) = 0.653131938 + 0.244542459 i + 0.366813688 j + 0.489084918 k

Note:

-In general, the hypergeometric series converge if  | a | < 1
 

4°)  Harmonic Numbers
 

 "HRMA" calculates the generalized harmonic numbers defined as  H_m(a) = 1 + 2 ^-a + 3 ^-a + ............ + m ^-a
 

Data Registers:           •  R00 = n > 1                ( Registers R00 thru Rnn are to be initialized before executing "HRMA" )

                                      •  R01   ......  •  Rnn = the n components of the anion

                                         Rn+1 ..........  R3n:  temp

  >>>  When the program stops:     R01   ......  Rnn = the n components of the result
 

Flags:  /
Subroutines:  "X^A"                                      ( cf "Anions for the HP41" )
                         "ST*A"  "DSA"          ( cf "Anionic Special Functions(I) for the HP41" )
 
 
 

 
    ( 81 bytes / SIZE 3n+1 )
 
 

   where  m is a positive integer and 1.nnn   is the control number of  the result

Example:    Calculate  H₁₂( 1 + 2 e₁ + 3 e₂ + 4 e₃ + 5 e₄ + 6 e₅ + 7 e₆ + 8 e₇ )                ( octonion ->  8  STO 00 )

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

    12  XEQ "HRMA"  >>>>  1.008                                ---Execution time = 44s---

   R01 = 0.248285998     R05 = 0.031438213
   R02 = 0.012575285     R06 = 0.037725855
   R03 = 0.018862927     R07 = 0.044013498
   R04 = 0.025150570     R08 = 0.050301140     whence

   H₁₂( 1 + 2 e₁ + 3 e₂ + 4 e₃ + 5 e₄ + 6 e₅ + 7 e₆ + 8 e₇ )  =  0.248285998 + 0.012575285 e₁ + 0.018862927 e₂ + 0.025150570 e₃
                                                                                           + 0.031438213 e₄ + 0.037725855 e₅ + 0.044013498 e₆ + 0.050301140 e₇
 

5°)  Bessel-Clifford Functions
 

Formula:    C_m(a) = SUM_{k=0,1,...,infinity}  ( 1/Gamma(k+m+1) )  a^k / k!
 

Data Registers:           •  R00 = n > 1                ( Registers R00 thru Rnn are to be initialized before executing "BSCLA" )

                                      •  R01   ......  •  Rnn = the n components of the anion

                                         Rn+1 ..........  R3n:  temp

  >>>  When the program stops:     R01   ......  Rnn = the n components of the result
 

Flags:  /
Subroutines:   "A*A1"                                               ( cf "Anions for the HP41" )
                         "ST*A"  "ST/A"  "DSA"          ( cf "Anionic Special Functions(I) for the HP41" )
                         "1/G+"                                          ( cf "Gamma function for the HP41" )
 
 

 
    ( 105 bytes / SIZE 3n+1 )
 
 

   where   1.nnn   is the control number of  the result

Example:     C_PI( 1 + 2 i + 3 j + 4 k ) = ?               ( quaternion ->  4  STO 00 )

     1   STO 01
     2   STO 02
     3   STO 03
     4   STO 04

    PI   XEQ "BSCLA"  >>>>    1.004                                           ---Execution time = 17s---

   R01 = 0.070702245
   R02 = 0.069832002
   R03 = 0.104748003
   R04 = 0.139664004

Whence,   C_PI( 1 + 2 i + 3 j + 4 k ) = 0.070702245 + 0.069832002 i + 0.104748003 j + 0.139664004 k

Notes:

-"BSCLA" does not work if m is a negative integer.
-But it will if you use the program HGFB listed in "Anionic Special Functions(II)"  paragraph 6-b),
-The routine may be simplified a lot:
 
 

6°)  Whittaker-M Functions
 

Formula:

    M_q,p(a) = exp(-a/2)  a^p+1/2  M( p-q+1/2 , 1+2p , a )               where  M = Kummer's functions
 

Data Registers:           •  R00 = n > 1                ( Registers R00 thru Rnn are to be initialized before executing "WHIMA" )

                                      •  R01   ......  •  Rnn = the n components of the anion

                                         Rn+1 ..........  R3n:  temp

  >>>  When the program stops:     R01   ......  Rnn = the n components of the result
 

Flags:  /
Subroutines:  "A*A1"  "A^X"   "E^A"              ( cf "Anions for the HP41" )
                         "ST/A"  "KUMA"          ( cf "Anionic Special Functions(I) for the HP41" )
 
 

 
    ( 91 bytes / SIZE 3n+1 )
 
 

   where   1.nnn   is the control number of  the result

Example:      q = sqrt(2)   p = sqrt(3)    a = 0.4 + 0.5 i + 0.6 j + 0.7 k               ( quaternion ->  4  STO 00 )

    0.4  STO 01
    0.5  STO 02
    0.6  STO 03
    0.7  STO 04

     2   SQRT
     3   SQRT   XEQ "WHIMA"  >>>>  1.004                                           ---Execution time = 21s---

    R01 = -0.807065423
    R02 =  0.373202939
    R03 =  0.447843527
    R04 =  0.522484115

-So,  M_q,p(a) = -0.807065423 + 0.373202939 i + 0.447843527 j + 0.522484115 k
 

7°)  Whittaker-W Functions
 

Formula:

  W_q,p(a) = [ Gam(2p) / Gam(p-q+1/2) ]  M_q,-p(a) + [ Gam(-2p) / Gam(-p-q+1/2) ]  M_q,p(a)     assuming  2p is not an integer.
 

Data Registers:           •  R00 = n > 1                ( Registers R00 thru Rnn are to be initialized before executing "WHIWA" )

                                      •  R01   ......  •  Rnn = the n components of the anion

                                         Rn+1 ..........  R4n:  temp      R4n+1 = p  ,  R4n+2 = q

  >>>  When the program stops:     R01   ......  Rnn = the n components of the result
 

Flags:  /
Subroutines:  "WHIMA"                                   ( cf paragraph 6 above )
                         "ST*A"                     ( cf "Anionic Special Functions(I) for the HP41" )
                         "1/G+"                                ( cf "Gamma function for the HP41" )
 
 

 
    ( 182 bytes / SIZE 4n+3 )
 
 

   where   2.p  is not an integer  &  1.nnn   is the control number of  the result

Example:      q = sqrt(2)    p = sqrt(3)    a = 0.4 + 0.5 i + 0.6 j + 0.7 k               ( quaternion ->  4  STO 00 )

    0.4  STO 01
    0.5  STO 02
    0.6  STO 03
    0.7  STO 04

     2   SQRT
     3   SQRT   XEQ "WHIWA"  >>>>  1.004                                           ---Execution time = 56s---

    R01 =  2.275449624
    R02 = -1.129250840
    R03 = -1.355101007
    R04 = -1.580951175

Whence,        W_q,p(a)  = 2.275449624 - 1.129250840 i - 1.355101007 j - 1.580951175 k

Note:

-"WHIWA" does not work if 2.p is an integer.
 

8°)  Generalized Error-Functions
 

"GERFA" calculates  erf_m (a) = a  exp ( -a^m )  M ( 1 ; 1+1/m ; a^m )          where   M = Kummer's function
 
 

Data Registers:           •  R00 = n > 1                ( Registers R00 thru Rnn are to be initialized before executing "GERFA" )

                                      •  R01   ......  •  Rnn = the n components of the anion

                                         Rn+1 ..........  R4n:  temp

  >>>  When the program stops:     R01   ......  Rnn = the n components of the result
 

Flags:  /
Subroutines:  "KUMA"   "ST*A"   "AMOVE"  "ASWAP"        ( cf "Anionic Special Functions(I) for the HP41" )
                         "A*A1"   "A^X"   "E^A"                                          ( cf "Anions for the HP41" )

-The M-Code routines may be replaced by the corresponding focal programs...
 
 

 
( 57 bytes / SIZE 4n+1 )
 
 

   Where  1.nnn   is the control number of  the result

Example:      m = sqrt(2)    a = 1.1 + 1.2 i + 1.3 j + 1.4 k               ( quaternion ->  4  STO 00 )

    1.1   STO 01
    1.2   STO 02
    1.3   STO 03
    1.4   STO 04

     2  SQRT   XEQ "GERFA"  >>>>   1.004                                           ---Execution time = 37s---

    R01 =  1.205349648
    R02 = -0.208378333
    R03 = -0.225743195
    R04 = -0.243108056

 Whence    erf_m(a) = 1.205349648 - 0.208378333 i - 0.225743195 j - 0.243108056 k

Notes:

-The anion a must remains "small"
-If m = 2, we get "the" error-function divided by  2/sqrt(pi)
-To get the "normalized" generalized error-function, multiply the result by  Gamma(m+1) / sqrt(PI)
 

9°)  Lambert W-Function
 

-Here we use Newton's method to solve the equation   w exp(w) = a    where  a  is a given anion # -1
 

Data Registers:           •  R00 = n > 1                ( Registers R00 thru Rnn are to be initialized before executing "LBWA" )

                                      •  R01   ......  •  Rnn = the n components of the anion

                                         Rn+1 ..........  R4n:  temp

  >>>  When the program stops:     R01   ......  Rnn = the n components of the result
 

Flags:  /
Subroutines:  AMOVE  ASWAP   ST*A   DSA                        ( cf "Anionic Special Functions(I) for the HP41" )
                         LNA  E^A  A*A1  A-A  1/A                                             ( cf "Anions for the HP41" )

-The M-Code routines may of course be replaced by focal programs
 
 

 
   ( 88 bytes / SIZE 4n+1 )
 
 

   Where  1.nnn   is the control number of  the result

Example:       a = 1 + 2 i + 3 j + 4 k               ( quaternion ->  4  STO 00 )

     1   STO 01
     2   STO 02
     3   STO 03
     4   STO 04

     XEQ "LBWA"  >>>>   1.004                                           ---Execution time = 38s---

     R01 = 1.281459811
     R02 = 0.304051227
     R03 = 0.456076840
     R04 = 0.608102453

  Whence    W(a) = 1.281459811 + 0.304051227 i + 0.456076840 j + 0.608102453 k

Notes:

-The real parts of the successive approximations are displayed ( line 12 )
-This routine does not work if  a = -1 but w(-1) = -0.3181315052 + 1.337235701 i
 

10°) Continued Fractions
 

 "CFRA" calculates the continued fraction  f(a) =  b₀ + ( a₁/b₁+ ) ( a₂/b₂+ ) ................... ( a_n/b_n+ ) ......   by the modified Lentz's algorithm,

   assuming that all the anions have proportional imaginary parts.
 
 

Data Registers:           •  R00 = n > 1                ( Registers R00 thru Rnn are to be initialized before executing "CFRA" )

                                      •  R01   ......  •  Rnn = the n components of the anion

                                         Rn+1 ..........  R7n:  temp

  >>>  When the program stops:     R01   ......  Rnn = the n components of the result
 

Flags:  /
Subroutines:     "AJ"  "BJ"    ( see below )
                            AMOVE  ASWAP  ST*A   DSA                        ( cf "Anionic Special Functions(I) for the HP41" )
                            A*A1  A+A  "1/A"                                                           ( cf "Anions for the HP41" )

-The M-Code routines may be replaced by focal programs

-"AJ" must take the anion a in R01 thru Rnn ,  j  in synthetic register M
  and returns the anion  a_j  in R01 thru Rnn without disturbing R00 and R2n+1 thru R7n
-Same thing for "BJ"
 
 

 
    ( 193 bytes / SIZE 7n+1 )
 
 

   Where  1.nnn   is the control number of  the result

Example1:       a = 1 + 2 i + 3 j + 4 k               ( quaternion ->  4  STO 00 )

     1   STO 01
     2   STO 02
     3   STO 03
     4   STO 04

 with the continued fraction defined as   b₀ = 1 ,  a_j = 2 a + n  ,  b_j = a² + n²     ( j > 0 )
 
 

 
   XEQ "CFRA"  >>>>      1.004                                           ---Execution time = 95s---

   R01 =  1.036327819
   R02 = -0.143096562
   R03 = -0.214644843
   R04 = -0.286193124

Whence,   f( 1 + 2 i + 3 j + 4 k ) = 1.036327819 - 0.143096562 i - 0.214644843 j - 0.286193124 k

Notes:

-The HP-41 displays the real parts of the successive approximations.
-Since we use A*A1, "CFRA" will not work if the anions have non-proportional imaginary parts.
-If, however, only b₀ does not satisfy this property, calculate f(a) with b₀ = 0 and add the "true" b₀ at the end

Example2:                                        Error-Function

-The error function may be computed by an ascending series for small arguments.
-But there is also a continued fraction that is preferable for large arguments, provided  Re(a) > 0

   1 - erf(a) = (1/sqrt(PI))  exp (-a² )  (1/a+) (1/2/a+) (1/a+) (3/2/a+ ) .....................
 
 

 
-To calculate  erf ( 1.8 + 1.9 i + 2 j + 2.1 k )               ( quaternion ->  4  STO 00 )

     1.8   STO 01
     1.9   STO 02
      2     STO 03
     2.1   STO 04

     XEQ "ERFA2"  >>>>         1.004                                           ---Execution time = 2m33s---

    R01 = -533.4595229
    R02 =  434.2164660
    R03 =  457.0699647
    R04 =  479.9234626

-So,   erf ( 1.8 + 1.9 i + 2 j + 2.1 k ) = -533.4595229 + 434.2164660 i + 457.0699647 j + 479.9234626 k
 

Notes:

-This continued fraction may also be calculated more directly ( from right to left )
  if you specify the number of terms N to be taken into account:
 
 

 
    ( 72 bytes / SIZE 2n+1 )
 
 

   Where  N is the number of terms in the continued fraction  and  1.nnn   is the control number of  the result

Example:    Calculate  erf ( 1.8 + 1.9 i + 2 j + 2.1 k )               ( quaternion ->  4  STO 00 )

     1.8   STO 01
     1.9   STO 02
      2     STO 03
     2.1   STO 04

-With   N = 19     XEQ "ERFA3"  >>>>         1.004                                           ---Execution time = 33s---

    R01 = -533.4595240
    R02 =  434.2164654
    R03 =  457.0699637
    R04 =  479.9234616

-So,   erf ( 1.8 + 1.9 i + 2 j + 2.1 k ) = -533.4595240 + 434.2164654 i + 457.0699637 j + 479.9234616 k

Notes:

-With  N = 20 , we get the same result.

-Of course, this is not the recommended way to evaluate a continued fraction,
 but roundoff-errors are probably smaller and moreover, SIZE 2n+1 is enough instead of SIZE 7n+1 with "CFRA"
-So, "ERFA3" may be used with 128-ons !
 

References:

[1]   Abramowitz and Stegun - "Handbook of Mathematical Functions" -  Dover Publications -  ISBN  0-486-61272-4
[2]   http://functions.wolfram.com/