The Museum of HP Calculators

# Error Function for the HP-41

This program is Copyright © 2006-2007 by Jean-Marc Baillard and is used here by permission.

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/pi1/2)  §0x  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/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

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/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°) Using 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 )

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