hp41programs

Error

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"
 02  2
 03  PI
 04  .147
 05  *
 06  /

 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/