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<Y? 32 GTO 01 |
33 31
34 AMOVE 35 RCL 00 36 E3 37 / 38 ISG X 39 CLD 40 END |
( 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<Y? 85 GTO 10 86 31 87 AMOVE 88 RCL 00 89 E3 90 / 91 ISG X 92 CLA 93 CLD 94 END |
( 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:
[1] Abramowitz and Stegun - "Handbook of Mathematical Functions"
- Dover Publications - ISBN 0-486-61272-4
[2] http://functions.wolfram.com/