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

( 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/