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

[Gil]

New version 7

- with some improvements
- and inclusion of other related functions.

1) Same input checking as in version 6.

2) Improvements on output "outlook".

3) Inverse functions, named now IERF and IERFc with capital I
(and no more iERF and iERFc , as small i has another meaning —> see 4) below).

4) Inclusion of 3 functions related to the error function:
- iERFc
- ERFi
- w.FADDEEVA-KRAMP
For definitions, see Wikipedia/French and Wolfram Alpha, or 10) for the codes.

5) Inverse_ERF(p)
Work with positive values of p only
IF p-1+1=p
Then Inverse_ERF(p) = ERFc(1-p)
(we want to be sure that we don't lose digits by calculating 1-p;
see, for instance, the results for p=5.0000000005E-4:
by Wolfram Alpha —> 4.4311349177240266E-4
by integral as p-1+1≠p —> 4.43113491773E-4
by built-in UTPN (in HP48-50G), but not using Inverse_ERFc(1-p) —> 4.43113492614E-4
by Inverse_ERFc(1-p) —> 4.43113491729E-4)
Else use Solver with initial value of p with always the Integral calculation, that gives better results.
END
IF initially p was negative
Then result =: -result
Else result = last, positive result
END

6) Inverse_ERFc(p)
IF p > 1
Then p:= 2-p
(we calculate initially for resulting positive x values)
Else p= initial value p
END
IF p now > 0.9 (it's an approximate value)
Then Solver with always the Integral calculation, that gives better results
(example results for x =.9999:
by Wolfram Alpha: 8.8622692777289E-5
by integral: 8.86226927774E-5
by built-in UTPN in HP48-50G: 8.86226927139E-5)
Else use the built-in UTPN function in HP48-50G.
END
IF you proceded initially to the change of p into 2-p
Then final result = -result
Else final result = same, positive result
END

7) ERF¦c, together with ERF (a help program used internally in ERF¦c, but not to be used directly)

Both programs still should be improved (some final digits induly may not appear in the outputs ERF¦c

8) All functions take one single argument.

9) All functions admit complex arguments on the form (a b), less the inverse functions IERF and IERFc (that take only real arguments).

10) Codes listing

ERF¦c
1 input: real x or a complex (a b)
2 outputs : the two corresponding "1-p" and "p" values ERFc and ERF
\\<< \"
1 Arg

.real 
x 
or
.complex 
(a b)
\" DROP DTAG EVAL 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
\\<<
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 .01 \\<= 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 ERFc 1 SWAP - *
END im 0 ==
IF
THEN ROT ERFc UNROT
ELSE ROT DROP DUP 1 SWAP - UNROT
END
END
\\>> 'IERR' PURGE
\\>>

IERF
1 real input
1 real output
\\<< \"
1 Arg

.a real 
a
, |
a
| \\<= 1
(
not
 a complex !)\" DROP DUP DUP IM 0 \\=/ SWAP RE ABS 1 > OR
IF
THEN DROP \"
Not
 OK as:
.|Real| > 1
or
.Complex with
IM\\=/0 or
|RE| > 1\" DOERR
END DUP ABS 1 >
IF
THEN \"|a|
must be \\<= 1\" DOERR
END DUP \"
Inv
ERF
(\" SWAP + \")\" + SWAP DUP 0 \\>= 1 -1 IFTE SWAP ABS \\-> p
\\<<
CASE p 1 - 1 + p ==
THEN 1 p - IERFc SWAP DROP
END p '\\.S(-x,x,e^-X^2.,X)/\\v/\\pi' - 'x' 2 ROOT { x IERR } PURGE
END *
\\>>
\\>>

IERFc
1 complex or real input
1 complex output
\\<< \"
1 Arg

.a real 
a
 in [0 2]
(
not
 a complex !)\" DROP \\-> a
\\<< \"
Inv
ERF
c

\" a TYPE 1 \\=/
IF
THEN \"(\" + a + \")\"
ELSE a
END + a IM 0 \\=/ a RE 0 < OR a RE 2 > OR
IF
THEN DROP \"
Not
 OK as:
.Real <0, >2
or
.Complex with
IM\\=/0 or
RE <0, >2
\" DOERR
ELSE PUSH -105 CF a 1 \\<=
IF
THEN 1
ELSE 2 a - 'a' STO -1
END
CASE a 0 ==
THEN '\\oo' EVAL
END a 1 ==
THEN 0
END a .9 \\<=
THEN
\\<< a 0 .5 x UTPN 2 * -
\\>> 'x' 2. ROOT
END 1 a - RE '\\.S(-x,x,e^-X^2.,X)/\\v/\\pi' - 'x' 2 ROOT
END * { x IERR } PURGE
END POP
\\>>
\\>>


iERFc
1 complex or real input
1 complex output
\\<< DUP \\-> z
\\<< \"
i
ERF
c

\" z TYPE 1 ==
IF
THEN z
ELSE \"(\" + z + \")\"
END + SWAP ERF\\166c 4 ROLL 3 DROPN z * z SQ NEG EXP \\pi \\v/ / \\->NUM - DUP NEG \"Wiki_French\" \\->TAG SWAP \"Wolfram\" \\->TAG
\\>>
\\>>

ERFi
1 input complex or real input
1 complex output
\\<< DUP \\-> z
\\<< \"ERF
i

\" z TYPE 1 ==
IF
THEN z
ELSE \"(\" + z + \")\"
END + SWAP i * \\->NUM ERF\\166c 4 ROLLD 3 DROPN i / \\->NUM
\\>>
\\>>

W.FADD
1 complex or real input
1 complex output
\\<< DUP \\-> z
\\<< \"wFADDEEVA-KRAMP
\" z TYPE 1 ==
IF
THEN z
ELSE \"(\" + z + \")\"
END + SWAP i NEG * \\->NUM ERF\\166c 4 ROLL 3 DROPN z SQ NEG EXP * \\->NUM
\\>>
\\>>

ERFc
Not for the user
Program used internally by ERF¦c
\\<< \"
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
\\>>

11) Request for Inverse_ERF (IERF) and Inverse_ERFc (IERFc)

For both above mentioned functions and without having to use an integral (which is quite time consuming if using the real calculator), I would duly appreciate to have corresponding, possible algorithms and codes allowing to get — in all possible situations — the same or better results.