The Museum of HP Calculators

Error function and it's inverse for the HP-32sII

This program is by Karl Bédard and is used here by permission.

This program is supplied without representation or warranty of any kind. Karl Bédard 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

There are 2 programs included, the error function and its inverse. The value of erf(z) can be found in mathematical tables, as build-in functions in calculators and spread sheets. The following approximation is from Professor Nathan Cheung, U.C. Berkeley:

erf(z)=1 - (a1×T + a2×T2 + a3×T3 + a4×T4 + a5×T5)×e(-z2)
where T=1/(1 + Pz) and other numerical values are:

P = 0.3275911,
a1 = 0.254829592,
a2 = - 0.284496736,
a3 = 1.421413741,
a4 = -1.453152027,
a5 =1.061405429.

THE FUNCTION ERF PROGRAM:

Program giving the numerical value of erf(z) when z is in the X stack. To call this function, simply press XEQ R. Should be accurate up to 1 part in 106. Therefore, not all displayed digits are significant. Use with care. You need the labels R, U and Q. This new version is more stable, precise and also works with the new hp33s. 
Example:
0.19
XEQ R (running...)
answer = 0.211839…

Main Label: 
01 LBL R
02 x≠0?
03 GTO U
04 0
05 RTN
CK D4FD 7.5 bytes
(similarly in HP33s, CK=059D LN=27)

because T is not behaving properly when z is negative: 
01 LBL U
02 x>0?
03 GTO Q
04 +/-
05 XEQ Q
06 +/-
07  RTN
CK AD1D 10.5 bytes
(similarly in HP33s, CK=00E0 LN=21)

Polynomial using Horner: 
01 LBL Q
02 STO Z
03 0.3275911
04 ×
05 1
06 +
07 1/X
08 STO T
09 1.061405429
10 ×
11 1.453152027
12 -
13 RCL×T
14 1.421413741
15 +
16 RCL×T
17 0.284496736
18 -
19 RCL×T
20 0.254829592
21 +
22 RCL×T
23 RCL Z
24 x2
25 +/-
26 ex
27 ×
28 +/-
29 1
30 +
31 RTN
CK C48E   94.5 bytes
(similarly in HP33s, CK=AB42 LN=189)

THE INVERSE OF THE FUNCTION ERF PROGRAM:

The inverse of that function can also be made from the following labels. This will allow you to get the value of z from erf(z). You will need the labels Z,J,T,V,C,O and R,U,Q from above. Should be accurate up to 1 part in 106. Therefore, not all displayed digits are significant. 
Example:
1
XEQ Z (running...)
answer = DIVIDE BY ZERO (what we want to see)

Also, 0.9999 will return something like 2,75. 
Let’s not forget that erf(z) is close enough to 1 when z=2,4.

Main Label: 
01 LBL Z
02 STO I
03 1
04 STO C
05 STO D
06 STO i
06 0
07 STO N
08 STO O
10 GTO C
11 RTN
CK AE92 16.5 bytes
(similarly in HP33s, CK=865E LN=57)

The inverse of erf requires the use of a program for linear
intrapolation. From points (STO C,STO D) and (STO N,STO O),
will give value of J for a point on the same line 
(STO I, J?). It comes as a bonus with the error fonction... 

01 LBL J
02 RCL I
03 RCL -C
04 RCL O
05 RCL-D
06 ×
07 RCL N
08 RCL-C
09 ÷
10 RCL+D
11 RTN
CK 8012  16.5 bytes
(similarly in HP33s, CK=DF86 LN=33)

Output 
01 LBL O
02 XEQ J
03 RTN
CK 796B 4.5 bytes
(similarly in HP33s, CK=0C10 LN=9)

01 LBL A
02 STO O
03 XEQ R
04 STO N
05 22
06 STO i
07 GTO C
08 RTN
CK 3637 12 bytes
(similarly in HP33s, CK=6E35 LN=36)

01 LBL V
02 STO D
03 XEQ R
04 STO C
05 1
06 STO i
07 GTO C
08 RTN
CK 2453 12 bytes
(similarly in HP33s, CK=7131 LN=36)

Convergence test: (squared is for negative values of N) 
01 LBL C
02 RCL I
03 RCL -N
04 x2
05 1E-12
06 x>y?
07 GTO O
08 XEQ J
09 GTO (i)
10 RTN
CK BBFD 23 bytes
(similarly in HP33s, CK=CB53 LN=42)
The convergence should take under 10 seconds.

Go back to the software library
Go back to the main exhibit hall