RE: HP48-HP50G : Error functions ERF(z), FRESNEL integrals, PROBIT, etc.

[Gil]

Version 9f

In the previous versions 9, the simplification of the input argument was incorrect when in presence of a complex arguments z=a+ib, with b≠0 (but everything was correctly calculated).

Now the 'arg.num' subroutine was duly fixed up.

Please note for INVERSE_ERF, called 'INV¦',
when using complex argument z=a+ib, b≠0,
that the results are quite accurate for abs(argument z) < 0.91.

Please note that this last uploaded version 9f
has a "print" error at variable VERSIO.

Please change it

from
Version 9e
2023.12.14
campart
@hotmail.com

into
Version 9f
2023.12.14
campart
@hotmail.com

[Gil]

Please note that the last uploaded version 9f
has a "print" error at variable VERSIO.

Please change it

from
Version 9e
2023.12.14
campart
@hotmail.com

into
Version 9f
2023.12.14
campart
@hotmail.com

ERF "\\<< \"
1 Arg

.real 
x

(no limit for real)
.complex
- 
(
a b
)

-
 '
a+i*b
'
- '
r*e^(i\\Gh)
'
- '
r*EXP(i\\Gh)
'

(b\\=/0 &

|
a
| 
& 
|
b
| 
\\<=9)
\" DROP DTAG DUP DUP DUP DUP PUSH euler DUP DUP RE SWAP IM \\oo \\->NUM 9 .01 \\-> \\<-z0 \\<-z \\<-re \\<-im \\<-max.x \\<-max.z critic
\\<< arg.num inf.erf max.arg DROP -105 SF EVAL ABS \\oo ==
IF
THEN SIGN -105 CF '\\oo' EVAL *
END DUP DUP \"ERF

(
\" SWAP DROP \\<-z0 + \"
)
\" + \"'\" \"\" SREPL DROP DUP \"F\" \"F
c
\" SREPL DROP 4 ROLLD SWAP DUP TYPE 9 ==
IF
THEN DROP \\<-z
END
CASE \\<-re DUP 0 \\>= 1 -1 IFTE SWAP ABS '\\oo' EVAL ==
THEN 1 * DUP 1 SWAP - ROT 5 ROLL DROP2 UNROT
END DROP \\<-im 0 == \\<-re ABS critic \\<= AND
IF
THEN DROP '\\.S(-\\<-re,\\<-re,e^-X^2,X)/\\v/\\pi' \\->NUM
ELSE \\<-im 0 == \\<-re 0 \\<= AND
IF
THEN -1
ELSE 1
END SWAP \\<-im 0 ==
IF
THEN ABS
END complex
END \\<-im 0 ==
IF
THEN \\<-re ABS critic \\<=
IF
THEN 1 SWAP
ELSE
END 4 ROLL DROP DUP 1 SWAP - ROT * DUP 0 <
IF
THEN DUP NEG 1 + ROT DROP SWAP
END \\<-re ABS critic >
IF
THEN SWAP
END UNROT
ELSE ROT DROP SWAP UNROT
END
END \\<-im 0 \\=/
IF
THEN 4 ROLL DROP
END
\\>> 'IERR' PURGE POP
\\>>" "



" INV¦ "\\<< \"
1 Arg

.a real 
p
, |
p
| \\<= 1
<or equivalent
complex 
z 
= 
p

- 
(
a b
)

-
 '
a+i*b
'
- '
r*e^(i\\Gh)
'
- '
r*EXP(i\\Gh)
' 
>


\" DROP PUSH DUP DUP 0 .9065 \\-> p0 \\<-z0 p \\<-z critic
\\<< p0 euler DUP DUP '\\<-z' STO 'p' STO arg.num \\<-z0 'p0' STO IM 0 \\=/ p ABS 1 \\>= AND p RE ABS 1 > OR
IF
THEN POP p0 \"
Not
 OK as:
.Real 
|
p
|
 > 1
or
.Complex p
with 
|
p
|
 > 1\" DOERR
END -105 CF \"
Inv
ERF
\" \"
(
\" + p0 + \"
)
\" + \"'\" \"\" SREPL DROP p IM 0 \\=/
IF
THEN p inv\\166z
ELSE p RE DUP 'p' STO 0 \\>= 1 -1 IFTE p ABS 'p' STO
CASE p 1 ==
THEN '\\oo' EVAL
END p 0 ==
THEN 0
END 1 p - critic \\<=
THEN
\\<< 1 p - 0 .5 x UTPN 2 * -
\\>> 'x' 2 ROOT
END p '\\.S(-x,x,e^-X^2,X)/\\v/\\pi' - 'x' 2 ROOT
END * { x IERR } PURGE
END POP
\\>>
\\>>" "



" INV¦c "\\<< \"
1 Arg

.a real 
p
 in [0 2]
<or equivalent
complex 
z 
= 
p

- 
(
a b
)

-
 '
a+i*b
'
- '
r*e^(i\\Gh)
'


 - '
r*EXP(i\\Gh)
' 
>
\" DROP PUSH DUP DUP 0 .9065 \\-> p0 \\<-z0 p \\<-z critic
\\<< p0 euler DUP '\\<-z' STO 'p' STO arg.num \\<-z0 'p0' STO \"
Inv
ERF
c

\" \"
(
\" + p0 + \"
)
\" + \"'\" \"\" SREPL DROP p IM 0 \\=/ 1 p - ABS 1 \\>= AND p RE 0 < OR p RE 2 > OR
IF
THEN DROP POP p0 \"
Not
 OK as:
.Real p<0, >2
or
.Complex with
IM(p)\\=/0 &

|
1-p
|
 \\>= 1
\" DOERR
ELSE 1 p - IM 0 \\=/
IF
THEN 1 p - inv\\166z
ELSE -105 CF p RE DUP 'p' STO 1 \\<=
IF
THEN 1
ELSE 2 p - 'p' STO -1
END
CASE p 0 ==
THEN '\\oo' EVAL
END p 1 ==
THEN 0
END p critic \\<=
THEN
\\<< p 0 .5 x UTPN 2 * -
\\>> 'x' 2 ROOT
END 1 p - RE '\\.S(-x,x,e^-X^2,X)/\\v/\\pi' - 'x' 2 ROOT
END * { x IERR } PURGE
END
END POP
\\>>
\\>>" "



" PROBIT "\\<< \"
1 Arg

.a real 
p
 in [0 1]
<or equivalent
complex 
z 
= 
p

- 
(
a b
)

-
 '
a+i*b
'
- '
r*e^(i\\Gh)
'
- '
r*EXP(i\\Gh)
' 
>


\" DROP PUSH DUP DUP 0 \\-> p0 \\<-z0 p \\<-z
\\<< p0 euler DUP DUP '\\<-z' STO 'p' STO arg.num \\<-z0 'p0' STO \"Probit
\" \"
(
\" + p0 + \"
)
\" + \"'\" \"\" SREPL DROP p IM 0 \\=/ p RE 1 > OR p RE 0 < OR
IF
THEN DROP POP \"
Not
 OK as:
.Real p<0, >1
or
.Complex with
IM(p) \\=/ 0 or
RE(p) <0, >1
\" DOERR
END SWAP DROP -105 CF
CASE p 1 ==
THEN \\oo EVAL
END p 0 ==
THEN \\oo EVAL NEG
END p RE 2 * 1 - INV\\166 2 \\v/ * \\->NUM SWAP DROP
END POP
\\>>
\\>>" "



" ERFi "\\<< \"
1 Arg

.real 
x


|
x
|
 \\<=9
.Or complex
- 
(
a b
)

-
 '
a+i*b
'
- '
r*e^(i\\Gh)
'
- '
r*EXP(i\\Gh)
'

 
|
a
| 
& 
|
b
| 
\\<=7
\" DROP PUSH DUP DUP DUP DUP euler DUP DUP RE SWAP IM 9 7 0 \\-> \\<-z0 \\<-z \\<-re \\<-im \\<-max.x \\<-max.z \\<-inf1x
\\<< arg.num inf2x DUP \"OK\" SAME
IF
THEN DROP inf1x.erfi \\<-inf1x 1 ==
IF
THEN inf1x.end \"ERF
i

\" SWAP + \"(+\" \"(\" SREPL 5 ROLL 5 ROLL 3 DROPN SWAP
ELSE max.arg -105 SF EVAL ABS \\oo ==
IF
THEN SIGN -105 CF '\\oo' EVAL * '\\<-z0' STO
ELSE DROP
END \"ERF
i

\" \"
(
\" + \\<-z0 + \"
)
\" + \"'\" \"\" SREPL ROT DROP2
CASE \\<-z ABS '\\oo' EVAL ==
THEN \\<-z SIGN -105 CF \\oo EVAL *
END \\<-z i * \\->NUM ERF 4 ROLLD 3 DROPN i / \\->NUM
END DUP IM 0 == { RE } IFT POP
END
ELSE \"ERF
i

\" SWAP + SWAP 5 ROLL 5 ROLL 5 ROLL 3 DROPN
END POP
\\>>
\\>>" "



" iERFc "\\<< \"
1 Arg

.real 
x

(no limit)
.or complex
- 
(
a b
)

-
 '
a+i*b
'
- '
r*e^(i\\Gh)
'
- '
r*EXP(i\\Gh)
'

(b\\=/0 &

|
a
| 
or 
|
b
| 
\\<=9)
\" DROP PUSH DUP DUP DUP DUP euler DUP DUP RE SWAP IM \\oo \\->NUM 9 0 \\-> \\<-z0 \\<-z \\<-re \\<-im \\<-max.x \\<-max.z \\<-inf1x
\\<< arg.num max.arg \\->NUM ABS \\oo \\->NUM ==
IF
THEN SIGN -105 CF '\\oo' EVAL * '\\<-z0' STO
ELSE DROP
END \"
i
ERF
c

\" \"
(
\" + \\<-z0 + \"
)
\" + \"'\" \"\" SREPL ROT DROP2
CASE \\<-z \\->NUM '\\oo' EVAL \\->NUM ==
THEN \"Undef by Wolfram
as '0/0', or \\->0
(try 
i
ERF
c
(x) with
x=10 x=20 x=33)\"
END \\<-z \\->NUM '\\oo' EVAL NEG \\->NUM ==
IF
THEN -105 CF \\oo EVAL NEG DUP NEG
ELSE \\<-z ERF ROT 4 ROLLD 3 DROPN \\<-z * \\<-z SQ NEG EXP \\pi \\v/ / - \\->NUM DUP IM 0 == { RE } IFT DUP NEG
END DUP 0 ==
IF
THEN DROP
ELSE \"Wiki_French\" \\->TAG SWAP \"Wolfram\" \\->TAG
END
END POP
\\>>
\\>>" "



" ERFcx "\\<< \"
1 Arg

.real 
x 


|
x
|
 \\<=33 or
.complex
- 
(
a b
)

-
 '
a+i*b
'
- '
r*e^(i\\Gh)
'
- '
r*EXP(i\\Gh)
'

(b\\=/0 &

|
a
| 
or 
|
b
| 
\\<=9)
\" DROP PUSH DUP DUP DUP DUP euler DUP DUP RE SWAP IM 33 9 0 \\-> \\<-z0 \\<-z \\<-re \\<-im \\<-max.x \\<-max.z \\<-inf1x
\\<< arg.num max.arg -105 SF EVAL ABS \\oo ==
IF
THEN SIGN -105 CF '\\oo' EVAL * '\\<-z0' STO
ELSE DROP
END \"ERF
cx

\" \"
(
\" + \\<-z0 + \"
)
\" + \"'\" \"\" SREPL ROT DROP2
CASE \\<-z '\\oo' EVAL ==
THEN 0
END \\<-z '\\oo' EVAL NEG ==
THEN '\\oo' EVAL
END \\<-z ERF ROT 4 ROLLD 3 DROPN \\<-z SQ EXP \\->NUM * DUP IM 0 == { RE } IFT
END POP
\\>>
\\>>" "



" FADDEEVA "\\<< \"
1 Arg

.real 
x


|
x
|
 \\<=9 or
.complex
- 
(
a b
)

-
 '
a+i*b
'
- '
r*e^(i\\Gh)
'
- '
r*EXP(i\\Gh)
'

(b\\=/0 &

|
a
| 
or 
|
b
| 
\\<=9)
\" DROP PUSH DUP DUP DUP DUP euler DUP DUP RE SWAP IM 9 9 0 \\-> \\<-z0 \\<-z \\<-re \\<-im \\<-max.x \\<-max.z \\<-inf1x
\\<< arg.num max.arg -105 SF EVAL ABS \\oo ==
IF
THEN SIGN -105 CF '\\oo' EVAL * '\\<-z0' STO
ELSE DROP
END \"
w
_
F
aaddeva-Kramp
\" \"
(
\" + \\<-z0 + \"
)
\" + \"'\" \"\" SREPL ROT DROP2 \\<-z i NEG * \\->NUM ERF ROT \\<-z SQ NEG EXP * \\->NUM 4 ROLLD 3 DROPN DUP IM 0 == { RE } IFT POP
\\>>
\\>>" "



" DAWSON "\\<< \"
1 Arg

.real 
x


|
x
|
 \\<=9 or
.complex
- 
(
a b
)

-
 '
a+i*b
'
- '
r*e^(i\\Gh)
'
- '
r*EXP(i\\Gh)
'
|
a
| 
& 
|
b
| 
\\<=9
 

\" DROP PUSH DUP DUP DUP DUP euler DUP DUP RE SWAP IM 9 9 0 \\-> \\<-z0 \\<-z \\<-re \\<-im \\<-max.x \\<-max.z \\<-inf1x
\\<< arg.num max.arg -105 SF EVAL ABS \\oo ==
IF
THEN SIGN -105 CF '\\oo' EVAL * '\\<-z0' STO
ELSE DROP
END \"
F
_Dawson
\" \"
(
\" + \\<-z0 + \"
)
\" + \"'\" \"\" SREPL ROT DROP2 \\<-z ABS \\->NUM '\\oo' EVAL \\->NUM ==
IF
THEN 0
ELSE \\<-z ERFi \\<-z SQ NEG EXP \\pi \\v/ * 2 / * \\->NUM SWAP DROP DUP IM 0 == { RE } IFT
END POP
\\>>
\\>>" "



" euler "\\<< -2. CF -103. SF -105. SF DUP \\->STR \"EXP\" \"e^\" SREPL DROP \\-> e s
\\<< s \"e^(\" POS DUP \\-> p
\\<< 0. ==
IF
THEN e
ELSE s p 2. + s SIZE SUB \"'\" SWAP + \"\\pi\" \"1\" SREPL DROP \"i\" \"1\" SREPL 0. ==
IF
THEN DROP e
ELSE OBJ\\-> \\->NUM DUP FP DUP 0. \\=/ SWAP .5 \\=/ AND
IF
THEN DROP e
ELSE \\->Q \\pi * -105. CF i * EXP EVAL p 2. \\=/
IF
THEN s 1. p 2. - SUB \"'\" + OBJ\\-> *
END
END
END
END \\->NUM DUP IM 0. ==
IF
THEN RE
END
\\>>
\\>>
\\>>" "



" arg.num "\\<< -105 SF \\<-z0 TYPE 9 == \\<-z0 TYPE 1 == \\<-z0 IM 0 == AND OR
IF
THEN \\<-z0 \"
=\" + \\<-z DUP IM 0 == \\<-z RE FP 0 == AND { R\\->I } IFT + '\\<-z0' STO
END
\\>>" "



" complex "\\<< DUP IM 0 ==
IF
THEN RE 0 .5 ROT UTPN 2 *
ELSE C\\->R SWAP DUP 0 ==
IF
THEN .00000000001 SWAP
ELSE 1
END \\-> y x z
\\<< '\\.S(-x,x,e^-X^2,X)/\\v/\\pi' \\->NUM 0 .5 x UTPN 2 * 'e^-x^2/(2*\\pi*x)*(1-COS(2*x*y)+i*SIN(2*x*y))' \\->NUM '2/\\pi*e^-x^2' \\->NUM '\\GS(k=1,30,e^(-k^2/4)/(k^2+4*x^2)*(2*x*(1-COS(2*x*y)*COSH(k*y))+k*SIN(2*x*y)*SINH(k*y)+i*(2*x*SIN(2*x*y)*COSH(k*y)+k*COS(2​*x*y)*SINH(k*y))))' \\->NUM * + NEG DUP 4 ROLL SWAP - UNROT + SWAP z 0 ==
IF
THEN C\\->R SWAP DROP 0 SWAP R\\->C SWAP C\\->R SWAP DROP 1 SWAP R\\->C SWAP
END
\\>>
END
\\>>" "



" max.arg "\\<< \\<-im 0. ==
IF
THEN \\<-re ABS \\<-max.x >
IF
THEN message
END
ELSE \\<-re ABS \\<-max.z > \\<-im ABS \\<-max.z > OR
IF
THEN message
END
END
\\>>" "



" message "\\<< DROP2 POP \\<-max.x \\oo \\->NUM ==
IF
THEN \"
Not
 OK with
z=a+ib, b\\=/0
if 
|
a
|

or 
|
b
| 
>\" \\<-max.z R\\->I + DOERR
ELSE \"
Not
 OK with
.Real
if 
|
x
| 
>\" \\<-max.x DUP \\oo \\->NUM ==
IF
THEN DROP \"\\oo\"
ELSE R\\->I
END + \"
\\183Or z=a+ib &
b\\=/0 if 
|
a
|

or 
|
b
| 
>\" + \\<-max.z R\\->I + DOERR
END
\\>>" "



" inf.erf "\\<< \\<-re ABS \\oo \\->NUM == \\<-im ABS \\oo \\->NUM == OR
IF
THEN DROP \\<-im ABS \\oo \\->NUM == \\<-re ABS \\oo \\->NUM == AND
IF
THEN \"Undef\" DUP \\<-re 0. \\>= \"(\\oo\" \"(-\\oo\" IFTE \\<-im 0. \\>= \"+i\\oo)\" \"-i\\oo)\" IFTE +
ELSE
CASE \\<-re ABS \\oo \\->NUM ==
THEN \\<-re 0. \\>= 1. -1. IFTE DUP 1. SWAP - SWAP \"(\" \\<-re 0. \\>= \"\" \"-\" IFTE + \"\\oo\" + \\<-im 0. \\=/
IF
THEN \\<-im 0. \\>= \"+i\" \"-i\" IFTE + \\<-im ABS +
END \")\" +
END \"1\" \\<-im 0. \\>= \"-\" \"+\" IFTE + \"i\\oo\" + DUP \"1\" \"\" SREPL DROP \"-\" \"\" SREPL 0. ==
IF
THEN \"+\" \"-\" SREPL DROP
END \"(\" \\<-re 0. \\=/
IF
THEN \\<-re + \" \" +
END \\<-im 0. \\>= \"+\" \"-\" IFTE + \"i\\oo)\" +
END
END \"ERF
\" SWAP + \"(+\" \"(\" SREPL DROP DUP \"F\" \"F
c
\" SREPL DROP SWAP ROT 4. ROLL UNROT 6. ROLL 6. ROLL DROP2 POP KILL
END
\\>>" "



" inf2x "\\<< \\<-im ABS \\oo \\->NUM == \\<-re ABS \\oo \\->NUM == AND
IF
THEN \"Undef\" \\<-re 0. \\>= \"(\\oo\" \"(-\\oo\" IFTE \\<-im 0. \\>= \"+i\\oo)\" \"-i\\oo)\" IFTE +
ELSE \"OK\"
END
\\>>" "



" inf1x.erfi "\\<< \\<-re ABS \\oo \\->NUM == \\<-im ABS \\oo \\->NUM == OR
IF
THEN DROP 1. '\\<-inf1x' STO \\<-im ABS \\oo \\->NUM == \\<-re ABS \\oo \\->NUM == AND
IF
THEN \"Undef\" \\<-im 0. \\>= \"+i\\oo)\" \"-i\\oo)\" IFTE
ELSE
CASE \\<-im \\oo \\->NUM ==
THEN \"i\" \"+i\\oo)\"
END \\<-im NEG \\oo \\->NUM ==
THEN \"-i\" \"-i\\oo)\"
END \\<-re 0. \\>= \"\\oo\" \"-\\oo\" IFTE \\<-im 0. \\=/
IF
THEN \\<-im 0. > \"+i\" \"-i\" IFTE \\<-im ABS +
ELSE \"\"
END \")\" +
END
END
END
\\>>" "



" inf1x.end "\\<<
CASE \\<-re 0. ==
THEN \"(\"
END \\<-re \\oo \\->NUM ==
THEN \"(\\oo\"
END \\<-re NEG \\oo \\->NUM ==
THEN \"(-\\oo\"
END \"(\" \\<-re + \" \" +
END SWAP +
\\>>" "



" inv¦z "\\<< { 1 } 100 \"This last value can
be incr. \\->160 or +
if 
|
z
|
 \\-> 1
\" DROP \\-> z c max
\\<< 1 max
FOR k 0 0 k 1 -
FOR m c m 1 + GET c k 1 - m - 1 + GET * m 1 + m 2 * 1 + * / + EVAL
NEXT c SWAP \\->NUM + 'c' STO
NEXT 0 0 max
FOR k c k 1 + GET 2 k * 1 + / \\pi \\v/ 2 / z * 2 k * 1 + ^ * \\->NUM +
NEXT
\\>>
\\>>" "



" EXPL "\\<< \"ERF(z) ERFi(z): 
OK

\\-> z=(a b),(\\177\\oo b),
 (a \\177 \\oo),(\\177\\oo \\177\\oo) 
iERFc(z) ERFcx
FADDEEVA(z) DAWSON(z):
\\-> z 
still \\=/
 \\177\\oo(i) 
INV\\166 INV\\166c PROBIT: 
OK

To speed up \\-> change

max
 in 'inv\\166z' to 
60
\" 1 DISP 7 FREEZE
\\>>" "



" VERSIO "Version 9f
2023.12.14

campart
@hotmail.com"

[Gil]

Version 10d
(correction of versions 10b & 10c, which were deleted;
see also note at the end of this post)

(version 10 non "published")


Thanks Albert Chan's insight, here is a simplified version of inverse_ERF (as usual called INV¦) & inverse_ERFc (as usual called INV¦c), using Newton's method in the subroutine 'inv¦z' when the argument p is a complex number z=a+ib, with b≠0.

However, for argument abs(p) a real number in [0 1], in case of inverse_ERF (called INV¦),
& for argument (p) a real number in [0 2], in case of inverse_ERFc (called INV¦c),

the usual formulae were left.

INV¦ "\\<< \"
1 Arg

.a real 
p
, |
p
| \\<= 1
<or equivalent
complex 
z 
= 
p

- 
(
a b
)

-
 '
a+i*b
'
- '
r*e^(i\\Gh)
'
- '
r*EXP(i\\Gh)
' 
>


\" DROP PUSH DUP DUP 0 .9065 \\-> p0 \\<-z0 p \\<-z critic
\\<< p0 euler DUP DUP '\\<-z' STO 'p' STO arg.num \\<-z0 'p0' STO IM 0 \\=/ p RE ABS 8.8 > p IM ABS 8.8 > OR AND p IM 0 == p RE ABS 1 > AND OR
IF
THEN POP p0 \"
Not
 OK as:
.Real 
|
p
|
 > 1
.or Complex
z=a+ib, b\\=/0
with 
|
a
|

or 
|
b
| 
>8.8\" DOERR
END -105 CF \"
Inv
ERF
\" \"
(
\" + p0 + \"
)
\" + \"'\" \"\" SREPL DROP p IM 0 \\=/
IF
THEN p DUP DUP inv\\166z
ELSE p RE DUP 'p' STO 0 \\>= 1 -1 IFTE p ABS 'p' STO
CASE p 1 ==
THEN '\\oo' EVAL
END p 0 ==
THEN 0
END 1 p - critic \\<=
THEN
\\<< 1 p - 0 .5 x UTPN 2 * -
\\>> 'x' 2 ROOT
END p '\\.S(-x,x,e^-X^2,X)/\\v/\\pi' - 'x' 2 ROOT
END * { x IERR } PURGE
END POP
\\>>
\\>>"


INV¦c "\\<< \"
1 Arg

.a real 
p
 in [0 2]
<or equivalent
complex 
z 
= 
p

- 
(
a b
)

-
 '
a+i*b
'
- '
r*e^(i\\Gh)
'


 - '
r*EXP(i\\Gh)
' 
>
\" DROP PUSH DUP DUP 0 .9065 \\-> p0 \\<-z0 p \\<-z critic
\\<< p0 euler DUP '\\<-z' STO 'p' STO arg.num \\<-z0 'p0' STO \"
Inv
ERF
c

\" \"
(
\" + p0 + \"
)
\" + \"'\" \"\" SREPL DROP p IM 0 \\=/ p RE ABS 8 > p IM ABS 8 > OR AND p IM 0 == p RE 0 < AND OR p IM 0 == p RE 2 > AND OR
IF
THEN DROP POP p0 \"
Not
 OK as:
.Real p<0, >2
.or Complex
z=a+ib, b\\=/0
with 
|
a
|

or 
|
b
|
 >8
\" DOERR
ELSE 1 p - IM 0 \\=/
IF
THEN p 1 p - DUP inv\\166z
ELSE -105 CF p RE DUP 'p' STO 1 \\<=
IF
THEN 1
ELSE 2 p - 'p' STO -1
END
CASE p 0 ==
THEN '\\oo' EVAL
END p 1 ==
THEN 0
END p critic \\<=
THEN
\\<< p 0 .5 x UTPN 2 * -
\\>> 'x' 2 ROOT
END 1 p - RE '\\.S(-x,x,e^-X^2,X)/\\v/\\pi' - 'x' 2 ROOT
END * { x IERR } PURGE
END
END POP
\\>>
\\>>"
" inv¦z "\\<<
DO \\-> p z
\\<< z DUP ERF 4 ROLLD 3 DROPN p - DUP \\-> f
\\<< '2/\\v/\\pi*EXP(-z^2)+f*z' / - \\->NUM p SWAP z
\\>>
\\>>
UNTIL OVER - DUP RE ABS .000000001 < SWAP IM ABS .000000001 < AND
END UNROT DROP2
\\>>"

Please note that, for precision purpose, erf(z=a+ib, with b≠0) is imited to the following argument condition: abs(a)<= 9 and abs (b)<= 9.
As the calculation of 'INV¦' use erf(z=a+ib), INV¦(z=a+ib, with b≠0) should be theorically limited to the following argument condition: abs(a)<= 9 and abs (b)<= 9.
But, in order to avoid an intermediary result > 9 with Newton's method (in subroutine 'inv¦z'), the present version 10c sets for INV¦(z=a+ib, with b≠0) the following more restrictive argument condition: abs(a)<= 8.8 and abs (b)<= 8.8.

The same applies to INV¦(z=a+ib, with b≠0), whose argument condition is now limited, in this version 10d, to:
abs(a)<= 8 and abs (b)<= 8.

[Gil]

New version 11b

All the explanations relative to the argument were somewhat changed.


Subroutine inv¦z was completed (in bold) for special cases :

inv¦z
\\<<
DO \\-> p z
\\<< z DUP ERF 4 ROLLD 3 DROPN p - DUP \\-> f
\\<< '2/\\v/\\pi*EXP(-z^2)+f*z' / - \\->NUM DUP DUP RE ABS 9 > SWAP IM ABS 9 > OR
IF
THEN ROT DROP \"
Failure
 as
Interm.Newton
result (
a
,
b
)=
\" SWAP 1 RND + \", &
|
a
| or |
b
| \\>=9 not OK in 
ERF
\" + POP DOERR \"To force a sol, try,
for max.z, 
99
 instead
of 
9
, 
here
 & in 
ERF
\" DROP
END p SWAP z
\\>>
\\>>
UNTIL OVER - DUP RE ABS .0000000001 < SWAP IM ABS .0000000001 < AND
END UNROT DROP2
\\>>

[Gil]

New version 11c

Slight changes (uniformization) in the DOERR messages.

Subroutine inv¦z was completed (in bold) for special cases.

The rest remains mostly the same for the functions
ERF & ERFc, their inverse INV¦ & INV¦c
& PROBIT ERFi iERFc ERFcx FADDEEVA DAWSON according to Wolfram Alpha or Wikipedia definitions.

As previously, all of them take a single argument input and accept, as argument input, a real or complex number, according to the explanation given at the beginning of each function.

inv¦z
\<<
DO 9 \-> p z \<-max.z
\<< z DUP ERF 4 ROLLD 3 DROPN p - DUP \-> f
\<< '2/\v/\pi*EXP(-z^2)+f*z' / - \->NUM DUP DUP RE ABS \<-max.z > SWAP IM ABS \<-max.z > OR
IF
THEN ROT DROP "To get here a sol,
try & let 
\<-max.z 
=

99
, instead of 
9
,
in 
inv
\166
z
 & in 
ERF
" SWAP "
Failure
 as
Interm.Newton
result (
a
,
b
)=
" SWAP 1 RND + ", &
|
a
|" + \<-max.z 10 < " or " "or" IFTE + "|
b
| \>=" + \<-max.z R\->I + " not OK in 
ERF
" + UNROT SWAP ROT POP DOERR
END p SWAP z
\>>
\>>
UNTIL OVER - DUP RE ABS .0000000001 < SWAP IM ABS .0000000001 < AND
END UNROT DROP2
\>>

[Gil]

New version: 12c

The 7 main changes:

1) As all other functions, PROBIT function accepts now, as argument, a complex number.
Reason: PROBIT function is defined as k*Inverse_ERF, and the function Inverse_ERF does accept — theorically — a complex number as argument.

2) As PROBIT, the functions Inverse_ERF (called 'INV¦') and Inverse_ERFc (called 'INV¦c') use Newton's recursive algorithm when in presence of a complex number argument. Therefore, in order to prevent endless loops, a counter (with a maximum loop by default set at 90) in the subroutine 'inv¦z' has been introduced in the DO... UNTIL.

3) The second criterium to stop the DO... UNTIL in the subroutine 'inv¦z' is relative to finding a converging, acceptable solution. Now, this criterium was "facilitated", with the critical difference (to interrupt the DO... UNTIL routine) between two consecutive complex numbers found (a_n, b_n) and (a_n+1, b_n+1) passing from abs(a_n - a_n+1) & abs(b_n - b_n+1) = 10^-10 to 10^-9.

4) To avoid name confusion, the previously named function iERFc was "more correctly" renamed i¹ERFc, with a special mention at the beginning of the program, as the initial letter i does not stand here for the inverse function (inverse_ERFc). See Wikipedia or Wolfram Alpha for definition and signs, Wolfram definition reversing the terms.

Note: in the menu, you will see only the first characters i¹ERF (without the final letter c) of its full, real name i¹ERFc.
Reminder: the inverse_ERFc function is called here 'INV¦c' (and called INV¦ for the inverse_ERF function).

5) Similarly to the function i¹ERFc, a new function i²ERFc was created. See Wikipedia or Wolfram Alpha for definition and signs, Wolfram definition reversing the terms.

Note 1: in the menu, you will see only the first characters i²ERFc (without the final letter c) of its full, real name i²ERFc.
Note 2: the output values of i²ERFc function could not be checked with another program, so be careful when using that function.

6) Some parameters might be exceptionally changed (such as <—max.x, <—max.z <—loop.max or critic). Therefore, their default values are specially reminded, in the beginning of each program, in form of a string, should you change some of the corresponding values and want afterwards to get back their original values.

7) Special attention was given to have the same number of POP and PUSH commands in the different programs, in order to always finish/interrupt a program with 'ENVSTACK' containing, at the very end, an empty list { }.

PS
Two "cosmetical" corrections in version 12c
vs precious (now deleted) version 12b.

[Gil]

New version 13

The error function directory
includes now also Owen's T function.

[Gil]

New version 14c

The directory includes now,
besides the following functions
ERF(z, z a complex number) & ERFc(z) & their inverses INV¦(z) & INV¦c(z)
PROBIT(z)
i¹ERFc(z)
i²ERFc(z)
ERFi(z)
ERFcx(z)
FADDEEVA(z)
DAWSON(z)
OWEN(h & a, but both arguments only real numbers)
the Fresnel integrals
FresnelS & FresnelC

with the following restriction:


1 Arg

.Real 
x
, |
x
| less or equal 10
. Or complex (x,
0)
with |
x
| less or equal 10


. Or complex
z=
a
+i
b
, 
b
not equal 0
with |
a
| & |
b
|
 less or equal 5

Arg
 here normalized:
't^2*pi/2'
instead of 't^2'

—> lim x—>±infinity = 1/2
instead of 1/2*sqrt (pi) / sqrt(2)

Complex form accepted for all function
(with the exception of OWEN's T function) :
- 
(
a b
)

-
 '
a+i*b
'
- '
r*e^(i*pi)
'
- '
r*EXP(i*pi)
'


In this version 14c,
OWEN_T function was also improved regarding its 2 arguments (calculation).

Most functions are restricted in the argument range
(a DOERR message with the corresponding explanation will appear
for non valid argument),
but all of them accept the special argument ± infinity.

Names/labels of subroutines called with small letters
& functions names mostly with capitals.

Program specially designed for EMU48.

[LdBeth]

I have written ERF and iERF for real numbers in SysRPL using Sergei Winitzki's A handy approximation for the error function and its inverse which I believe the inverse erf function is as good as the UserRPL version posted here, and is much faster.

I have originally planned to use that to computing Q function and inverse Q function, but the accuracy is only as good as 3-4 digits after decimal point even with extended real been used, so I decided to give up until I figure out a cheap way to compute the function using Taylor series.

You may find the binary for HP50g (and I believed should be compatible with HP49g at least), with the source code

https://github.com/LdBeth/eflib/releases/tag/v0.9.0

[Gil]

New version 14 e

Now 2 ways of calculating the Fresnel integrals Fresnel_S & Fresnel_C:

a) as in previous version, the normalzed form,
with always 1 single input:
- the function argument x (real) or z (complex) in level 1;

b) new, since this version 14e, the Not normalized version,
with then 2 inputs:
- the letter N in stack level 2 (with N for Not normalized)
- the function argument x (real) or z (complex) in l stack evel 1.

[Gil]

New version 15 b

Now normalized Fresnel Integrals for real argument x
or z=(x, zero), with no limit for x.

[Gil]

New version 16

The subroutine with the calculation of FRESNEL Integrals for very large real numbers
('fresnelFASTxLARGE') contained severe errors that were duly fixed up in this version.

My apologises for the inconvenience.

Now the Fresnel results given by the HP50G are very similar with the ones computed by Wolfram.

Besides, in this new version, you can always compute both the normalized and Not normalized Fresnel Integrais F_C & F_S
(just adding, when calculating the Non normalized Fresnel Integrals, the letter N [N for Non normalized] in stack level 2, before the argument x or z always in stack level 1).

Regarding always the Fresnel integrals:
when argument x is a real number : any value from -infinity to +infinity;
when argument z is a complex number a+ib, b≠0: abs(a) & abs(b) less or equal to 5

[Gil]

New version 16e

Minor changes, mostly in the (result) presentation.

In particular, the subroutine 'fresnelFASTxLARGE' for the calculation of FRESNEL Integrals when dealing with very large real numbers
was renamed 'FRESNEL.xLARGE.fast', with capital letters at the beginning, as it can be used as an independent function (but only for real numbers x, with abs(x) in [4 infinity).

However, to avoid confusion, this function or subroutine 'FRESNEL.xLARGE.fast' was put at the last page of the ERROR directory, as the main function FRESNEL.z to be used normally (at page 2 of the error directory) will automatically access the above mentioned subroutine as soon as x > 10.

[Gil]

New version 16.08

3 not at all important changes:

a) clearer explanations at the beginning of
- FRESNEL.z main program
- & subroutine function FRESNEL...;

b) change of the name of the above subroutine into:
- 'fresnelFASTxLARGE' (at the end, preference was given to small initial letters, indicating clearly that fresnel is a subroutine);

c) of course, FRESNEL.z refers now consequently to
new name fresnelFASTxLARGE.

Though meaningless or useless, here are some results found by the HP50G programs that might be compared with Wolfram Alpha:

{"ERF
c


(
33
)
" 1.93206244519E-475
"ERF

(
33
)
" 1. "

" "
Inv
ERF
c


(
1.93206244519E-475
)
" 33. "

" "Probit

(
1
)
" '+infinity' "

" "
i¹
ERF
c


(
5
)
" :Wiki (Wolframopp): 1.4813429336E-13 "

" "
i¹
ERF
c


(
+infinity 
)
" "Undef by Wolfram
as '0/0' (or —>0:
try 
i¹
ERF
c
(x) with
x=10 x=20 x=33)" "

" "
i¹
ERF
c


(
-infinity 
)
" :Wiki (Wolfram—>opp): '+infinity ' "

" "
i²
ERF
c


(
25
)
" :Wiki (Wolfram—>opp): 3.3068803E-277 "

" "
normaliz
.FRESNEL
S


(
10
)
" .468169978589 "

normaliz
.FRESNEL
C


(
10
)
" .499898694222 "

" "fresnelFASTxLARGE" "
normaliz
.FRESNEL
S


(
5000*e^(pi*i)
=-5000
)
" -.499936338031 "


normaliz
.FRESNEL
C


(
5000*e^(pi*i)
=-5000
)
" -.500000002602 "

" "fresnelFASTxLARGE" "
normaliz
.FRESNEL
S


(
500000000000
)
" .5 "


normaliz
.FRESNEL
C


(
500000000000
)
" .500000000001 "

" "
normaliz
.FRESNEL
S


(
-infinity 
)
" '-1/2' "


normaliz
.FRESNEL
C


(
-infinity 
)
" '-1/2' "

" "
not
_norm.FRESNEL
S


(
+infinity 
)
" 'sqrt(2*pi)/4' "


not
_norm.FRESNEL
C


(
+infinity 
)
" 'sqrt(2*pi)/4' "

" "ERF
i


(
(4.,2.)
)
" (-20442.123626,3999.26730293) "

" "ERF
i


(
4+2*i
=(4.,2.)
)
" (-20442.123626,3999.26730293) "

" "ERF
cx


(
(5.,0.)
=5
)
" .110704637733 "

" "
w
_
F
aaddeva-Kramp

(
5
)
" (1.3887943865E-11,.11524596183) "

" "
F
_Dawson

(
5
)
" .102134074424 "

" "P
(
X>
h
 & 0<Y<
a
X
)
:
Owen_T

(
h=5,
a=3
)
" 1.4332578594E-7

Enjoy the use of the differents function with your HP50G phone application!