The Museum of HP Calculators

Error Function for the HP-41

This program is supplied without representation or warranty of any kind. Jean-Marc Baillard and The Museum of HP Calculators therefore assume no responsibility and shall have no liability, consequential or otherwise, of any kind arising from the use of this program material or any part thereof.

Overview

1°) Series Expansion & Continued Fraction ( slightly improved )
2°) A Series Expansion
3°) Using Kummer's Function

-This function is defined by: erf x = (2/pi^1/2) §₀^x e^{-t^2} dt

1°) Series Expansion & Continued Fraction

-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/pi^1/2) SUM_{n=0,1,2,.....} (-1)ⁿ x²ⁿ⁺¹ / (n! (2n+1))

2.e^{x^2} §_x^infinity 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

Example: 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/Pi^1/2) exp(-x²) SUM_{n=0,1,2,.....} ( 2ⁿ x²ⁿ⁺¹ )/[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°) Using Kummer's Function

-An even shorter routine may be written if you have loaded "KUM" in your HP-41

Formula: erf x = (2x/Pi^1/2) exp(-x²) M( 1 , 3/2 , x² )

Data Registers: R00 = x² , 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 )

Reference: Abramowitz and Stegun , "Handbook of Mathematical Functions" - Dover Publications - ISBN 0-486-61272-4

Go back to the HP-41 software library
Go back to the general software library
Go back to the main exhibit hall