# hp41programs

Anionic Functions3

# 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)! / am + m! / (2.am+1) + SUMk=1,2,.... B2k (2k+m-1)! / (2k)! / a2k+m       where  B2k are Bernoulli numbers

•   If m = 0  ,  y (a)  ~  Ln a - 1/(2a) - SUMk=1,2,.... B2k / (2k) / a2k              ( digamma function )

and the recurrence relation:   y(m) (a+p) = y(m) (a) + (-1)m m! [ 1/am+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" )

 01  LBL "PSINA"   02  STO M   03  RCL 00   04  .1   05  %   06  ISG X   07  +   08   E3   09  /   10  1   11  +   12  STO N   13  REGMOVE   14  RCL 00   15  3   16  *   17  .1   18  %   19  ISG X   20  +   21  RCL 00   22  -   23  CLRGX   24  LBL 01   25  RCL N   26  REGSWAP   27  REGMOVE   28  8   29  RCL M   30  +   31  RCL 01   32  X>Y?   33  GTO 00   34  1   35  LASTX   36  +   37  CHS   38  XEQ "A^X"   39  XEQ "DSA"   40  RCL 00   41  1 42  ST+ Y   43  ST+ IND Y   44  GTO 01   45  LBL 00   46  RCL 00   47   E3   48  /   49  ISG X   50  RCL 00   51  3   52  *   53  +   54   E3   55  /   56  1   57  +   58  REGMOVE   59  XEQ "A*A1"   60  XEQ "1/A"   61  RCL N   62  REGMOVE   63  RCL M   64  9   65  +   66  FACT   67  39.6   68  /   69  XEQ "ST*A"   70  XEQ "A*A1"   71  RCL M   72  7   73  +   74  FACT   75  ST- 01   76  XEQ "A*A1"   77  40   78  XEQ "ST/A"   79  RCL M   80  5   81  +   82  FACT 83  ST+ 01   84  XEQ "A*A1"   85  42   86  XEQ "ST/A"   87  RCL M   88  3   89  +   90  FACT   91  ST- 01   92  XEQ "A*A1"   93  60   94  XEQ "ST/A"   95  RCL M   96  1   97  +   98  FACT   99  ST+ 01 100  XEQ "A*A1" 101  6 102  XEQ "ST/A" 103  RCL 00 104  .1 105  % 106  ISG X 107  + 108  RCL 00 109 + 110   E3 111  / 112  1 113  + 114  REGSWAP 115  RCL M 116  FACT 117  XEQ "ST*A" 118  X<>Y 119  REGSWAP 120  RCL N 121  REGMOVE 122  RCL 00 123   E3 124  / 125  ISG X 126  LASTX 127  / 128  1 129  + 130  RCL 00 131  3 132  * 133  + 134  STO O 135  REGMOVE 136  XEQ "1/A" 137  RCL M 138  FACT 139  XEQ "ST*A" 140  RCL 00 141  ENTER^ 142  ST+ Y 143  LBL 03 144  RCL IND Y 145  RCL IND Y 146  + 147  2 148  / 149  STO IND Y 150  RDN 151  DSE Y 152  DSE X 153  GTO 03 154  RCL M 155  X=0? 156  GTO 00 157  1 158  - 159  FACT 160  ST+ 01 161  RCL N 162  REGMOVE 163  RCL O 164  REGMOVE 165  RCL M 166  CHS 167  XEQ "A^X" 168  XEQ "A*A1" 169  GTO 02 170  LBL 00 171  RCL N 172  REGMOVE 173  RCL O 174  REGMOVE 175  XEQ "LNA" 176  RCL 00 177  ENTER^ 178  ST+ Y 179  LBL 04 180  RCL IND Y 181  RCL IND Y 182  - 183  STO IND Y 184  RDN 185  DSE Y 186  DSE X 187  GTO 04 188  LBL 02 189  XEQ "DSA" 190  RCL O 191  RCL 00 192  - 193  REGMOVE 194  SIGN 195  CHS 196  RCL M 197  Y^X 198  CHS 199  XEQ "ST*A" 200  RCL 00 201   E3 202  / 203  ISG X 204  CLA 205  END

( 401 bytes / SIZE 4n+1 )

 STACK INPUT OUTPUT X m 1.nnn

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

Example1:    Calculate  y(3) ( 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(3) ( 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) = ( am / 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" )

 01  LBL "IGFA"  02  RCL 00  03  3  04  *  05  1  06  +  07  X<>Y  08  STO IND Y     09  1  10  STO Z 11  +  12  XEQ "KUMA"  13   E3  14  /  15  1  16  +  17  RCL 00  18  +  19  STO M  20  REGMOVE 21  RCL 00  22  3  23  *  24  1  25  +  26  RCL IND X     27  STO N  28  XEQ "A^X"  29  RCL M  30  REGSWAP 31  SIGN  32  CHS  33  XEQ "ST*A"  34  XEQ "E^A"  35  XEQ "A*A1"    36  RCL M  37  RCL 00  38  .1  39  %  40  + 41  +  42  REGMOVE  43  XEQ "A*A1"    44  RCL N  45  XEQ "ST/A"  46  X<>Y  47  CLA  48  END

( 103 bytes / SIZE 3n+2 )

 STACK INPUT OUTPUT X m 1.nnn

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:    Ba(p,q) = ( ap / 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" )

 01  LBL "IBFA"  02  CHS  03  1  04  +  05  RCL Y  06  LASTX  07  +  08  CF 04  09  XEQ "HGFA"  10  R^  11  STO M  12  X<>Y  13   E3  14  /  15  1  16  +  17  RCL 00  18  +  19  REGSWAP  20  X<>Y  21  XEQ "A^X"  22  XEQ "A*A1"  23  RCL M  24  XEQ "ST/A"  25  X<>Y  26  CLA  27  END

( 61 bytes / SIZE 3n+1 )

 STACK INPUTS OUTPUTS Y p p X q 1.nnn

where   1.nnn   is the control number of  the result

Example:    Compute  B0.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,    B0.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  Hm(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" )

 01  LBL "HRMA"  02  STO M  03  1  04  CHS  05  XEQ "ST*A"  06  RCL 00  07   E3 08  /  09  ISG X  10  STO N  11  LASTX          12  /  13  1  14  + 15  RCL 00  16  +  17  STO O  18  REGSWAP   19  RCL N  20  RCL 00  21  ST+ X 22  .1  23  %  24  +  25  +  26  CLRGX         27  LBL 01   28  RCL O 29  REGMOVE   30  RCL M  31  XEQ "X^A"  32  XEQ "DSA"    33  DSE M  34  GTO 01  35  RCL O 36  RCL 00  37  +  38  REGMOVE    39  RCL N  40  CLA  41  END

( 81 bytes / SIZE 3n+1 )

 STACK INPUT OUTPUT X m 1.nnn

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

Example:    Calculate  H12( 1 + 2 e1 + 3 e2 + 4 e3 + 5 e4 + 6 e5 + 7 e6 + 8 e7 )                ( 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

H12( 1 + 2 e1 + 3 e2 + 4 e3 + 5 e4 + 6 e5 + 7 e6 + 8 e7 )  =  0.248285998 + 0.012575285 e1 + 0.018862927 e2 + 0.025150570 e3
+ 0.031438213 e4 + 0.037725855 e5 + 0.044013498 e6 + 0.050301140 e7

5°)  Bessel-Clifford Functions

Formula:    Cm(a) = SUMk=0,1,...,infinity  ( 1/Gamma(k+m+1) )  ak / 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" )

 01  LBL "BSCLA"  02  CLA  03  STO M  04  RCL 00  05  .1  06  %  07  STO Z  08  ISG X  09  +  10   E3 11  /  12  1  13  +  14  REGMOVE   15  0  16  XEQ "ST*A"   17  SIGN  18  STO 01  19  RDN  20  + 21  STO N  22  REGMOVE   23  LBL 01  24  XEQ "A*A1"   25  RCL M  26  RCL O  27  1   28  +  29  STO O  30  ST+ Y 31  *  32  XEQ "ST/A"   33  XEQ "DSA"  34  X#0?  35  GTO 01  36  RCL N  37  REGSWAP     38  SIGN  39  RCL M  40  + 41  XEQ "1/G+"  42  XEQ "ST*A"   43  RCL 00  44   E3  45  /  46  ISG X  47  CLA  48  END

( 105 bytes / SIZE 3n+1 )

 STACK INPUT OUTPUT X m 1.nnn

where   1.nnn   is the control number of  the result

Example:     CPI( 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,   CPI( 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:

 01  LBL "BSCLA"  02  1  03  +  04   E-3  05  CHS  06  XEQ "HGFB"  07  END

6°)  Whittaker-M Functions

Formula:

Mq,p(a) = exp(-a/2)  ap+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" )

 01  LBL "WHIMA"  02  ST- Y  03  ST+ X  04  1  05  +  06  .5  07  RCL Z  08  - 09  X<>Y  10  XEQ "KUMA"  11   E3  12  /  13  1  14  +  15  RCL 00  16  + 17  STO M  18  REGMOVE      19  CLX  20  2  21  /  22  XEQ "A^X"  23  FRC  24  ST+ X 25  RCL M  26  +  27  REGSWAP      28  XEQ "A*A1"  29  RCL M  30  REGSWAP  31  RCL 00  32  + 33  REGMOVE      34  2  35  CHS  36  XEQ "ST/A"  37  XEQ "E^A"  38  XEQ "A*A1"  39  CLA  40  END

( 91 bytes / SIZE 3n+1 )

 STACK INPUTS OUTPUTS Y q / X p 1.nnn

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,  Mq,p(a) = -0.807065423 + 0.373202939 i + 0.447843527 j + 0.522484115 k

7°)  Whittaker-W Functions

Formula:

Wq,p(a) = [ Gam(2p) / Gam(p-q+1/2) ]  Mq,-p(a) + [ Gam(-2p) / Gam(-p-q+1/2) ]  Mq,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" )

 01  LBL "WHIWA"  02  RCL 00  03  4  04  *  05  2  06  +  07  X<> Z  08  STO IND Z  09  DSE Z  10  X<>Y  11  STO IND Z  12  XEQ "WHIMA"  13  RCL 00  14  4  15  *  16  2  17  +  18  RCL IND X  19  STO N 20  DSE Y  21  RCL IND Y  22  STO M  23  ST+ X  24  CHS  25  XEQ "1/G+"  26  STO O  27  .5  28  RCL M  29  RCL N  30  +  31  -  32  XEQ "1/G+"  33  RCL O  34  /  35  XEQ "ST*A"      36  RCL 00  37   E3  38  / 39  ISG X  40  LASTX  41  /  42  1  43  +  44  RCL 00  45  3  46  *  47  +  48  REGSWAP  49  RCL 00  50  -  51  REGMOVE  52  RCL N  53  RCL M  54  CHS  55  XEQ "WHIMA"  56  STO N  57  RCL 00 58  4  59  *  60  2  61  +  62  RCL IND X  63  CHS  64  DSE Y  65  RCL IND Y   66  STO M  67  +  68  .5  69  +  70  XEQ "1/G+"       71  X<> M  72  ST+ X  73  XEQ "1/G+"  74  RCL M  75  X<>Y  76  / 77  XEQ "ST*A"  78  4  79  RCL 00  80  ST* Y  81  LBL 01  82  RCL IND Y   83  RCL IND Y  84  +  85  STO IND Y        86  RDN  87  DSE Y  88  DSE X  89  GTO 01  90  X<> N  91  CLA  92  END

( 182 bytes / SIZE 4n+3 )

 STACK INPUTS OUTPUTS Y q / X p 1.nnn

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,        Wq,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  erfm (a) = a  exp ( -am )  M ( 1 ; 1+1/m ; am )          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...

 01  LBL "GERFA"  02  STO M  03  14  04  AMOVE  05  X<>Y  06  XEQ "A^X"  07  SIGN  08  RCL M  09  1/X  10  1  11  +  12  XEQ "KUMA"  13  12  14  ASWAP  15  SIGN  16  CHS  17  ST*A  18  XEQ "E^A"  19  A*A1  20  42  21  AMOVE  22  A*A1  23  END

( 57 bytes / SIZE 4n+1 )

 STACK INPUTS OUTPUTS X m 1.nnn

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    erfm(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

 01  LBL "LBWA"  02  14  03  AMOVE  04  SIGN  05  ST+ 01  06  XEQ "LNA"  07  13  08  AMOVE 09  LBL 01  10  31  11  AMOVE  12  VIEW 01  13  SIGN  14  CHS  15  ST*A  16  XEQ "E^A" 17  42  18  AMOVE         19  A*A1  20  32  21  AMOVE  22  A-A  23  12  24  ASWAP 25  SIGN  26  ST+ 01  27  XEQ "1/A"      28  A*A1  29  DSA  30   E-8  31  X

( 88 bytes / SIZE 4n+1 )

 STACK INPUTS OUTPUTS X / 1.nnn

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) =  b0 + ( a1/b1+ ) ( a2/b2+ ) ................... ( an/bn+ ) ......   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  aj  in R01 thru Rnn without disturbing R00 and R2n+1 thru R7n
-Same thing for "BJ"

 01  LBL "CFRA"  02  16  03  AMOVE  04  CLA  05  XEQ "BN"  06  XEQ 00  07  13  08  AMOVE  09  15  10  AMOVE  11  CLX  12  ST*A  13  14  14  AMOVE  15  GTO 10  16  LBL 00  17  XEQ 01  18  X#0?  19  RTN 20   E-41  21  STO 01  22  RTN  23  LBL 01  24  RCL 00  25  0  26  LBL 02  27  RCL IND Y  28  X^2  29  +  30  DSE Y  31  GTO 02  32  RTN  33  LBL 10  34  ISG M  35  CLX  36  61  37  AMOVE  38  XEQ "AN" 39  17  40  AMOVE  41  61  42  AMOVE  43  XEQ "BN"    44  14  45  ASWAP  46  72  47  AMOVE  48  A*A1  49  42  50  AMOVE  51  A+A  52  XEQ 00  53  XEQ "1/A"  54  14  55  AMOVE  56  51  57  AMOVE 58  XEQ "1/A"    59  27  60  ASWAP  61  A*A1  62  72  63  AMOVE  64  A+A  65  XEQ 00  66  15  67  AMOVE  68  42  69  AMOVE  70  A*A1  71  32  72  AMOVE  73  13  74  AMOVE  75  A*A1  76  VIEW 01 77  13  78  ASWAP        79  SIGN  80  ST- 01  81  XEQ 01  82  SQRT  83  2 E-9  84  X

( 193 bytes / SIZE 7n+1 )

 STACK INPUTS OUTPUTS X / 1.nnn

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   b0 = 1 ,  aj = 2 a + n  ,  bj = a2 + n2     ( j > 0 )

 01  LBL "AJ"  02  2  03  ST*A  04  RCL M  05  ST+ 01  06  RTN  07  LBL "BJ"  08  RCL M  09  X=0?  10  GTO 00  11  12  12  AMOVE  13  A*A1  14  RCL M  15  X^2  16  ST+ 01  17  RTN  18  LBL 00  19  ST*A  20  SIGN  21  STO 01  22  RTN  23  END

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 b0 does not satisfy this property, calculate f(a) with b0 = 0 and add the "true" b0 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 (-a2 )  (1/a+) (1/2/a+) (1/a+) (3/2/a+ ) .....................

 01  LBL "ERFA2"  02  XEQ "CFRA"  03  61  04  AMOVE  05  12  06  AMOVE  07  A*A1  08  SIGN  09  CHS 10  ST*A  11  XEQ "E^A"     12  32  13  AMOVE  14  A*A1  15  PI  16  SQRT  17  CHS  18  ST/A 19  SIGN  20  ST- 01  21  X<>Y  22  RTN  23  LBL "AJ"         24  CLX  25  ST*A  26  RCL M  27  1 28  X#Y?  29  GTO 00  30  STO 01  31  RTN  32  LBL 00           33  -  34  2  35  /  36  STO 01 37  RTN  38  LBL "BJ"         39  RCL M  40  X#0?  41  RTN  42  ST*A  43  RTN  44  END

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

 01  LBL "ERFA3"  02  STO M  03  12  04  AMOVE  05  LBL 01  06  XEQ "1/A"  07  RCL M 08  1  09  -  10  STO M  11  X=0?  12  GTO 00         13  2  14  / 15  ST*A  16  A+A  17  GTO 01  18  LBL 00  19  12  20  ASWAP         21  2 22  XEQ "A^X"    23  SIGN  24  CHS  25  ST*A  26  XEQ "E^A"  27  A*A1  28  PI 29  SQRT  30  CHS  31  ST/A  32  SIGN  33  ST- 01           34  X<>Y  35  END

( 72 bytes / SIZE 2n+1 )

 STACK INPUTS OUTPUTS X N 1.nnn

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:

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