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

[Gil]

I changed the way of calculating erf,
now with the corresponding integral (much longer, of course), instead of using the UTPN command.

Reason :
Try 0. 5 1E-9 UTPN —> .499999999436
× 2 —> .999999998872
erf=1-.999999998872 =.000000001128.
Some digits got logically lost by this procedure.

Now, with integral from -x to x gives the expected answer in our case : 1.12837916709E-9.

When executing
.000000001 ERF¦c

We get two results:

"ERFc
(.000000001)"
.999999998872

"ERF
(.000000001)"
1.12837916709E-9

No more necessary to use ERFc, but of course you should not delete it, as it is used in ERF¦C.

The input can be a real or a COMPLEX number.

Note that iERFc will, logically, not give back the initial value x

ERF¦c
\\<< \"
1 Arg

.real 
x 
or
.complex 
(a b)
\" DROP DUP DUP DUP TYPE 1. ==
IF
THEN \"ERF
\" SWAP +
ELSE \"ERF
(\" SWAP + \")\" +
END DUP \"F\" \"Fc\" SREPL DROP 4. ROLLD SWAP DUP DUP IM SWAP RE \\-> im re
\\<< im 0. == re ABS .01 \\<= AND
IF
THEN DROP 0. re 'e^-X^2.' X \\.S 2. * \\pi \\v/ / \\->NUM
ELSE im 0. == re 0. \\<= AND
IF
THEN -1.
ELSE 1.
END SWAP im 0. ==
IF
THEN ABS
END ERFc 1. SWAP - *
END im 0. ==
IF
THEN ROT ERFc UNROT
ELSE ROT DROP DUP 1. SWAP - UNROT
END
\\>> 'IERR' PURGE
\\>>

iERF
\\<< \"
1 Arg

.a real 
a
, |
a
| < 1
(
not
 a complex !)\" DROP DUP IM 0. ==
IF
THEN RE
ELSE \"Complex with IM\\=/0:
NOT allowed for INV\" DOERR
END DUP ABS 1. \\>=
IF
THEN \"|a|
must be < 1\" DOERR
END DUP \"iERF
(\" SWAP + \")\" + SWAP '\\.S(-A,A,e^-X^2.,X)/\\v/\\pi' - 'A' 2. ROOT { A IERR } PURGE
\\>>

iERFc
\\<< \"
1 Arg

.a real 
a
, 
a
 < 2
(
not
 a complex !)\" DROP DUP IM 0. ==
IF
THEN RE
ELSE \"Complex with IM\\=/0:
NOT allowed for INV\" DOERR
END DUP ABS 2. \\>=
IF
THEN \"|a|
must be < 2\" DOERR
END DUP \"iERF
(\" SWAP + \")\" + SWAP 1. - NEG '\\.S(-A,A,e^-X^2.,X)/\\v/\\pi' - 'A' 2. ROOT { A IERR } PURGE
\\>>

ERFc
\\<< \"
1 Arg

.real 
x 
or
.complex 
(a b)
\" DROP 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
\\<< 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*C​OS(2.*x*y)*SINH(k*y))))' \\->NUM * + - z 0. ==
IF
THEN C\\->R SWAP DROP 1. SWAP R\\->C
END
\\>>
END
\\>>