Error Function for the HP-41
Overview
1°) Series Expansion & Continued Fraction
2°) A Series Expansion
3°) Via Kummer's Function
4°) Asymptotic Expansion
5°) Inverse Error
Function
a) Ascending Series
b) Via Halley's
Method
6°) Approximate Formulae
7°) Generalized Error
Functions
1°) Series Expansion & Continued Fraction
-This function is defined by: erf x = (2/pi1/2) §0x e-t^2 dt
-To achieve the best accuracy,
the following program calculates erf x by a series expansion
for x < Ln 6
and 1 - erf x by a continued fraction otherwise (
x non-negative )
Formulae: erf x = (2/pi1/2) SUMn=0,1,2,..... (-1)n x2n+1 / (n! (2n+1))
2.ex^2 §xinfinity e-t^2
dt = 1/(x+0.5/(x+1/(x+1.5/(x+2/(x+....)))))
Data Registers: /
Flags: /
Subroutines: /
-Synthetic registers M , N , O may of course be replaced by any
data registers.
| 01 LBL "ERF"
02 ABS 03 SIGN 04 LASTX 05 ENTER^ 06 STO M 07 STO O 08 6 09 LN 10 X<=Y? 11 GTO 02 12 RDN 13 X^2 14 STO N 15 LBL 01 16 CLX |
17 RCL N
18 RCL O 19 * 20 CHS 21 R^ 22 ISG T 23 CLX 24 / 25 STO O 26 R^ 27 ST+ X 28 DSE X 29 / 30 X<>Y 31 ST+ Y 32 X#Y? |
33 GTO 01
34 ST+ X 35 PI 36 SQRT 37 / 38 ENTER^ 39 CHS 40 1 41 + 42 X<>Y 43 GTO 04 44 LBL 02 45 CLX 46 24 47 STO N 48 R^ |
49 +
50 X<>Y 51 ST+ X 52 / 53 LBL 03 54 + 55 RCL N 56 X<>Y 57 ST+ X 58 / 59 DSE N 60 GTO 03 61 + 62 1/X 63 X<>Y 64 X^2 |
65 CHS
66 E^X 67 * 68 PI 69 SQRT 70 / 71 ENTER^ 72 CHS 73 1 74 + 75 LBL 04 76 RCL M 77 SIGN 78 RDN 79 CLA 80 END |
( 110 bytes / SIZE 000 )
| STACK | INPUTS | OUTPUTS |
| Y | / | 1-erf x |
| X | x >= 0 | erf x |
| L | / | x |
Examples:
0.9 XEQ
ERF >>>> erf(0.9) = 0.7969082124
X<>Y 1-erf(0.9) = 0.2030917876
( in 7.6 seconds )
2.7
R/S >>>> erf(2.7)
= 0.9998656673 X<>Y
1-erf(2.7) = 0.0001343327399
( in 7.7 seconds )
2°) A Series Expansion
-If you don't need to know 1 - erf x very accurately for large
arguments, this second program is shorter ( but slower ) than "ERF"
-This routine uses the formula: erf x = (2/Pi1/2)
exp(-x2) SUMn=0,1,2,..... ( 2n x2n+1
)/[1.3.5.....(2n+1)]
Data Registers: /
Flags: /
Subroutines: /
| 01 LBL "ERF2"
02 1 03 RCL Y 04 ENTER^ 05 X^2 06 STO M 07 LBL 01 |
08 X<> T
09 RCL M 10 ST+ X 11 * 12 R^ 13 SIGN 14 R^ |
15 +
16 STO T 17 ST+ L 18 X<> L 19 / 20 STO T 21 X<>Y |
22 ST+ Y
23 X#Y? 24 GTO 01 25 ST+ X 26 0 27 X<> M 28 E^X |
29 PI
30 SQRT 31 * 32 / 33 END |
( 55 bytes / SIZE 000 )
| STACK | INPUTS | OUTPUTS |
| X | x | erf x |
Examples:
1 XEQ "ERF2" >>>> erf 1 = 0.842700793
( in 6 seconds )
2 R/S
>>>> erf 2 = 0.995322265 ( 11s
)
3 R/S
>>>> erf 3 = 0.999977910 ( 15s
)
4 R/S
>>>> erf 4 = 0.999999985 ( 22s
)
3°) Via Kummer's Function
-An even shorter routine may be written if you have loaded "KUM" in your HP-41
Formula: erf x = (2x/Pi1/2)
exp(-x2) M( 1 , 3/2 , x2 )
Data Registers:
R00 = x2 , R01 = 1 , R02 = 3/2 , R03 = x
Flags: /
Subroutine: "KUM" ( cf "Kummer's
Functions for the HP-41" )
| 01 LBL "ERF3"
02 STO 03 03 1.5 04 STO 02 05 SIGN 06 STO 01 07 X<>Y 08 X^2 09 XEQ "KUM" 10 RCL 03 11 ST+ X 12 * 13 RCL 00 14 E^X 15 PI 16 SQRT 17 * 18 / 19 END |
( 35 bytes / SIZE 004 )
| STACK | INPUTS | OUTPUTS |
| X | x | erf x |
Example:
2 XEQ "ERF3" >>>> erf 2 = 0.995322265 ( in 14s )
Note:
-The routine may be improved further with the M-Code routine "HGF+"
( cf "Hypergeometric Functions" )
| 01 LBL "ERF3"
02 STO 00 03 1.5 04 STO 02 05 SIGN 06 ENTER^ 07 STO 01 08 RCL 00 09 X^2 10 HGF+ 11 LASTX 12 CHS 13 E^X 14 * 15 RCL 00 16 ST+ X 17 * 18 PI 19 SQRT 20 / 21 END |
( 34 bytes / SIZE 003 )
| STACK | INPUTS | OUTPUTS |
| X | x | erf x |
Example:
2 XEQ "ERF3" >>>> erf 2 = 0.995322265
---Execution time = 6s---
4°) Asymptotic Expansion
-This program computes the complementary error function:
erfc x = 1 - erf x ~ exp(-x2)
/ (x sqrt(PI)) [ 1 + SUM i=1,2,..... (-1/(2x2))k
1.3.....(2k-1) ] for large positive x
Data Registers:
R00 to R03: temp
Flags: /
Subroutines: /
| 01 LBL "AERFC"
02 ABS 03 STO 00 04 SIGN 05 STO 01 06 STO 02 07 CLX 08 STO 03 |
09 LBL 01
10 RCL 02 11 CHS 12 RCL 03 13 1 14 + 15 STO 03 16 ST+ X |
17 DSE X
18 * 19 RCL 00 20 X^2 21 ST+ X 22 / 23 STO 02 24 RCL 01 |
25 +
26 STO 01 27 LASTX 28 X#Y? 29 GTO 01 30 RCL 00 31 PI 32 SQRT |
33 *
34 / 35 RCL 00 36 X^2 37 CHS 38 E^X 39 * 40 END |
( 54 bytes / SIZE 004 )
| STACK | INPUTS | OUTPUTS |
| X | x > 0 | erfc x |
where erfc x = 1 - erf x
Example:
4.91 XEQ "AERFC" >>>> 3.817216231 10 -12 ---Execution time = 14s---
Notes:
-The series diverges too soon for x = 4.9
-As usual, the M-Code routine HGF+ may be used for faster results.
-The asymptotic expansion is re-written:
erfc x = 1 - erf x ~ exp(-x2)
/ (x sqrt(PI)) 2F0 ( 1 , 1/2 ;; -1/x2
)
| 01 LBL "AERFC"
02 ABS 03 STO 00 04 1 05 STO 01 |
06 .5
07 STO 02 08 1/X 09 0 10 RCL 00 |
11 X^2
12 CHS 13 STO 03 14 1/X 15 HGF+ |
16 RCL 03
17 E^X 18 * 19 RCL 00 20 PI |
21 SQRT
22 * 23 / 24 END |
( 36 bytes / SIZE 004 )
| STACK | INPUTS | OUTPUTS |
| X | x > 0 | erfc x |
where erfc x = 1 - erf x
Example:
4.91 XEQ "AERFC" >>>> 3.817216229
10 -12
---Execution time = 7s---
5°) Inverse Error Function
a) Ascending Series
Formula:
erf -1(x) = SUM k=0 to infinity [ ck/(2k+1) ] [ x sqrt(PI) / 2 ]2k+1
with c0 = 1 and ck
= SUM m=0 to k-1 cm ck-1-m / (m+1)
/ (2m+1)
Data Registers:
R00 = c0 = 1 , R01 = c1 , ............... , Rkk =
ck , ...................
Flags: /
Subroutines: /
| 01 LBL "IERF"
02 CLA 03 X^2 04 PI 05 * 06 4 07 / 08 STO N 09 SQRT 10 STO M 11 1 12 STO 00 |
13 X<>Y
14 LBL 01 15 RCL O 16 E3 17 / 18 STO P 19 CLX 20 LBL 02 21 RCL O 22 RCL P 23 INT 24 - |
25 RDN
26 RCL IND T 27 RCL IND L 28 * 29 RCL P 30 INT 31 ISG X 32 CLX 33 ST/ Y 34 ST+ X 35 DSE X 36 / |
37 +
38 ISG P 39 GTO 02 40 STO IND P 41 ISG O 42 CLX 43 RCL O 44 ST+ X 45 ISG X 46 CLX 47 / 48 RCL N |
49 ST* M
50 CLX 51 RCL M 52 * 53 + 54 X#Y? 55 GTO 01 56 CLA 57 END |
( 91 bytes / SIZE maximum )
| STACK | INPUTS | OUTPUTS |
| X | x | erf -1 x |
Examples:
0.4 XEQ "IERF" >>>> erf
-1 ( 0.4 ) = 0.370807159
---Execution time = 55s---
0.7
R/S >>>> erf
-1 ( 0.7 ) = 0.732869078
---Execution time = 4m52s---
Notes:
-Obviously, this program is slow and uses many registers.
-SIZE 028 is required for the second example.
-If x is very close to 1, a root-finding program with "ERF" is by far
preferable:
b) Via Halley's Method
-This variant uses Halley's 3rd order method to solve f(x) = erf(x) - y = 0
xn+1 = xn + 2 f(xn) f '(xn) / [ 2 ( f '(xn) )2 - f(xn) f '(xn) ]
-The convergence is cubic.
-To achieve a better accuracy, we employ erf(x) - 1 when erf(x)
is close to 1
Data Registers:
R00 thru R05: temp ( R05 = x )
Flag: F01
Subroutine: "ERF" ( cf paragraph 1 )
| 01 LBL "IERF"
02 STO 05 03 ABS 04 STO 03 05 STO 04 06 CF 01 07 .99 08 X<Y? 09 SF 01 10 LBL 01 |
11 RCL 04
12 XEQ "ERF" 13 FS? 01 14 GTO 02 15 RCL 03 16 - 17 GTO 03 18 LBL 02 19 CLX 20 RCL 03 |
21 1
22 - 23 + 24 CHS 25 LBL 03 26 STO Y 27 RCL 04 28 ST* Y 29 X^2 30 CHS |
31 E^X
32 ST+ X 33 PI 34 SQRT 35 / 36 + 37 / 38 ST- 04 39 VIEW 04 40 ABS |
41 E-8
42 X<Y? 43 GTO 01 44 RCL 04 45 RCL 05 46 SIGN 47 * 48 CF 01 49 CLD 50 END |
( 79 bytes / SIZE 006 )
| STACK | INPUTS | OUTPUTS |
| Z | / | ~ 0 |
| Y | / | E-8 |
| X | x | erf -1 x |
With -1 < x < +1
Examples:
0.7
XEQ "IERF" >>>>
erf -1 ( 0.7 )
= 0.732869079
---Execution time = 16s---
0.999999
R/S >>>>
erf -1 ( 0.999999 ) = 3.458910737
---Execution time = 86s---
.9999999999 R/S
>>>> erf -1 ( 0.9999999999 ) = 4.572824967
---Execution time = 122s---
Notes:
-Line 39 VIEW 04 displays the absolute values of the successive
approximations.
-But the output has the correct sign !
-Z-output is erf ( erf -1 x ) - x and is therefore
- almost - zero.
6°) Approximate Formulae
-If a great accuracy is not needed, these short routines allow to compute
very fast the error function and its inverse.
Formulae: cf reference [3]
erf x ~ sgn(x) { 1 - exp [ -x2 ( 4/PI + a x2 ) / ( 1 + a x2 ) ] } where a = 0.147
erf -1 x ~ sgn(x) { -2/(a.PI) - (1/2) Ln ( 1
- x2 ) + [ - (1/a) Ln ( 1 - x2 ) + [ 2/(a.PI) + (1/2)
Ln ( 1 - x2 ) ]2 ]1/2 }1/2
where a = 0.147 too
Data Registers: /
Flags: /
Subroutines: /
| 01 LBL "ERF"
02 ENTER^ 03 X^2 04 .147 05 * 06 4 07 PI 08 / 09 X<>Y 10 ST+ Y 11 1 |
12 +
13 / 14 X<>Y 15 X^2 16 * 17 CHS 18 E^X 19 1 20 X<>Y 21 - 22 SQRT |
23 X<>Y
24 SIGN 25 * 26 END 01 LBL "IERF"
|
07 X<>Y
08 X^2 09 CHS 10 LN1+X 11 2 12 / 13 ST+ Y 14 .0735 15 / 16 CHS 17 RCL Y |
18 X^2
19 + 20 SQRT 21 X<>Y 22 - 23 SQRT 24 X<>Y 25 SIGN 26 * 27 END |
( 83 bytes / SIZE 000 )
| STACK | INPUTS | OUTPUTS |
| Y | / | x |
| X | x | erf x or erf -1 x |
Examples:
1.4 XEQ "ERF" >>>>
erf (1.4) = 0.9524
0.9 XEQ "IERF" >>>
erf -1 (0.9) = 1.1630
Notes:
-Execution time is about 1.5 second.
-With erf x , the relative precision of the formula is better
than 1.3 10 -4 for all reals x
-With erf -1 x , the relative precision of the formula is
better than 2 10 -3 for -1 < x < 1
7°) Generalized Error Functions
-The error function may be generalized by the formula:
erfn x = §0x exp(-tn) dt
and we could use a series expansion to compute this function,
at least for small arguments.
-But it may also be calculated via the incomplete gamma function by
erfn x = (1/n) igam((1/n) , xn)
-The following program employs the equivalent method: erfn
x = x exp(-xn) 1F1 ( 1 , 1+1/n ; xn
)
Data Registers: R00 thru R03: temp
Flags: /
Subroutine: "KUM" ( cf "Kummer's
Function for the HP-41" )
| 01 LBL "ERFN"
02 STO 03 03 X<>Y 04 Y^X 05 LASTX 06 1/X 07 1 08 STO 01 09 + 10 STO 02 11 X<>Y 12 XEQ "KUM" 13 RCL 00 14 CHS 15 E^X 16 * 17 RCL 03 18 * 19 END |
( 32 bytes / SIZE 004 )
| STACK | INPUTS | OUTPUTS |
| X | x | erfn x |
where erfn x = §0x exp(-tn) dt
Example:
2 SQRT
PI XEQ "ERFN" >>>> erfsqrt(2)
(PI) = 0.907310876
---Execution time = 17s---
Note:
-It will not be a surprise: here again, we can use HGF+ too !
Data Registers: R00 thru R02: temp
Flags: /
Subroutines: /
| 01 LBL "ERFN"
02 STO 00 03 X<>Y 04 Y^X 05 LASTX 06 1/X 07 1 08 STO 01 09 + 10 STO 02 11 LASTX 12 ENTER^ 13 R^ 14 HGF+ 15 LASTX 16 CHS 17 E^X 18 * 19 RCL00 20 * 21 END |
( 31 bytes / SIZE 003 )
| STACK | INPUTS | OUTPUTS |
| X | x | erfn x |
where erfn x = §0x exp(-tn) dt
Example:
2 SQRT
PI XEQ "ERFN" >>>> erfsqrt(2)
(PI) = 0.907310876
---Execution time = 8s---
Notes:
-You can add a normalization factor so that erf2 x = erf
x
-Thus, some authors define erfn x = §0x
exp(-tn) dt ( n! / sqrt(PI) )
-On the other hand, removing the normalization factor
for erf x is also possible to save a few bytes ...
References:
[1] Abramowitz and Stegun , "Handbook of Mathematical Functions"
- Dover Publications - ISBN 0-486-61272-4
[2] http://dlmf.nist.gov/7.16
[3] Sergei Winitzki - A handy approximation for the error
function and its inverse - February 6, 2008
[4] http://functions.wolfram.com/