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

[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.