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

[Gil]

Summary

Now, both ERF and ERFc functions allow as single argument
a real x (the returned result will be a real "probability")
or a complex number in the form (a b) (the returned result will be a complex number);

Their inverse functions iERF and iERFc, however, only admit here as single argument a real number/"probability" (the returned result will be the real number x).

Note
For real argument/numbers, the results calculated here for ERF and Erfc may differ by 1 unit in the last digit in comparison with the built-in functions erf(x) and erfc(x) of the HP PRIME.
For complex arguments, the results calculated here for ERF and ERFc may slightly differ in their last two last digits (the HP PRIME is then more accurate).

ERF
\\<< \"
1 Arg

.real 
x 
or
.complex 
(a b)
\" DROP DUP TYPE 1 \\=/
IF
THEN 0
ELSE C\\->R
END \\-> x y
\\<< x DUP SIGN SWAP ABS 0 .5 ROT UTPN -2 * 1 + * y 0 \\=/
IF
THEN '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 * + \\->NUM
END
\\>>
\\>>

iERF
\\<< \"
1 Arg

.a real
(
not
 a complex !)\" DROP 1 - -2 / NOVAL \\-> p X
\\<<
\\<< p 0 .5 X UTPN -
\\>> 'X' 2 # 2F12Ch SYSEVAL \"x of Erf\" \\->TAG
\\>>
\\>>
If you have another calculator than the HP50G, please don't use # 2F12Ch SYSEVAL.

Use instead:
iERF
\<< "
1 Arg

.a real
(
not
 a complex !)" DROP 1 - -2 / \-> p
\<<
\<< p 0 .5 X UTPN -
\>> 'X' 2 ROOT "x of Erf" \->TAG 'X' PURGE
\>>
\>>

ERFc
\\<< \"
1 Arg

.real 
x 
or
.complex 
(a b)
\" DROP ERF NEG 1 +
\\>>

iERFc
\\<< \"
1 Arg

.a real
(
not
 a complex !)\" DROP 2 / NOVAL \\-> p X
\\<<
\\<< p 0 .5 X UTPN -
\\>> 'X' 2 # 2F12Ch SYSEVAL \"x of Erfc\" \\->TAG
\\>>
\\>>
If you have another calculator than the HP50G, please don't use # 2F12Ch SYSEVAL.

Use instead:
iERFc
\<< "
1 Arg

.a real
(
not
 a complex !)" DROP 2 / \-> p
\<<
\<< p 0 .5 X UTPN -
\>> 'X' 2 ROOT "x of Erfc" \->TAG 'X' PURGE
\>>
\>>

[Gil]

Check the following:

ERF(+5.1) :.999999999999

ERF (+5.2'):1
(with the HP50G roundings; in fact,correct is exactly 0.99999999999980751)

'ERFc (up to now) : 0 (1-"1" above)
(in fact should be 1.9249061099972305*10^-13)

But of course iERFc(0) does not give the expected value of 5.2.

We have to try and get iERFc(ERFc(x)) = x (or a near value).

That was the point of Albert Chan to which I did not pay attention.

Finally, I changed the programs,
starting by calculating ERFc, then setting ERF=1-ERFc.

ERF
\\<< \"
1 Arg

.real 
x 
or
.complex 
(a b)
\" DROP ERFc NEG 1 +
\\>>

iERF
\\<< \"Only HP 50G


1 Arg

.a real
(
not
 a complex !)\" DROP 1 - -2 / NOVAL \\-> p X
\\<<
\\<< p 0 .5 X UTPN -
\\>> 'X' 2 # 2F12Ch SYSEVAL \"x of Erf\" \\->TAG
\\>>
\\>>

ERFc
\\<< \"
1 Arg

.real 
x 
or
.complex 
(a b)
\" DROP DUP IM 0 ==
IF
THEN RE 0 .5 ROT UTPN 2 *
ELSE C\\->R \\-> x y
\\<< 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 * + -
\\>>
END
\\>>

iERFc
\\<< \"Only HP 50G


1 Arg

.a real
(
not
 a complex !)\" DROP 2 / NOVAL \\-> p X
\\<<
\\<< p 0 .5 X UTPN -
\\>> 'X' 2 # 2F12Ch SYSEVAL \"x of Erfc\" \\->TAG
\\>>
\\>>

MiERFc
\\<< \"
1 Arg

.a real
(
not
 a complex !)\" DROP 2 / \\-> p
\\<<
\\<< p 0 .5 X UTPN -
\\>> 'X' 2 ROOT \"x of Erfc\" \\->TAG 'X' PURGE
\\>>
\\>>

MiERF
\\<< \"
1 Arg

.a real
(
not
 a complex !)\" DROP 1 - -2 / \\-> p
\\<<
\\<< p 0 .5 X UTPN -
\\>> 'X' 2 ROOT \"x of Erf\" \\->TAG 'X' PURGE
\\>>
\\>>


Now, with the new version 2:
ERFc(5.2) :1.92490611E-13
and iERFc (+1.92490611E-13) : 5.2 (exactly).

[Gil]

Slight improvement for

ERF
\\<< \"
1 Arg

.real 
x 
or
.complex 
(a b)
\" DROP DUP SIGN SWAP ABS ERFc NEG 1 + *
\\>>

[Albert Chan]

(11-09-2023 12:46 PM)Albert Chan Wrote:  Erfc(x) = « 0 0.5 ROT UTPN 2 * »

Getting Erfc from UTPN is accurate. (only scaling, no subtraction cancellation)
With good Erfc, we expected getting IErfc should also be accurate.

Unfortunately, we don't have same formula pattern for Erf(x)
(is there a normal distribution area function for center part?)

Here, we use erf(x) slope to improve Erf and IErf

erf(x)' = 2/sqrt(pi) * exp(-x*x)

Code:

function ierf_taylor(x)
    if x<0 then return -ierf_taylor(-x) end
    local y = ierfc(1-x)    -- assume good ierfc
    local h = x - (x-1+1)
    return h==0 and y or y + sqrt(pi)/2 * exp(y*y) * h
end
function erf_taylor(x)
    local y = 1 - erfc(x)   -- assume good erfc
    local h = x - ierf_taylor(y)
    return h==0 and y or y + 2/sqrt(pi) * exp(-x*x) * h
end

lua> x = 1e-9
lua> erf(x), 1-erfc(x)
1.1283791670955127e-09     1.1283791678806665e-09
lua> erf_taylor(x)
1.1283791670955127e-09

lua> ierf(x), ierfc(1-x)
8.862269254527581e-10       8.862269003885489e-10
lua> ierf_taylor(x)
8.86226925452758e-10

[Gil]

Again, a lesson...

Yes, I saw and appreciate the huge improvement for x very small.

Sorry for my ignorance, but I did not understand the code to try and adapt it to HP50G.

Unfortunately.

I try here
y = ierfc(1-x) -- assume good ierfc
With x=10^(-9) —> y=:x of Erfc: 8.86019642276E-10
local h = x - (x-1+1) —> should be zero?
return h==0 and y or y + sqrt(pi)/2 * exp(y*y) * h
But y small, y² very small & exp(y²) = 1
y=y+sqrt (pi)/2 *0 * 0 = y+0 = y
end
function erf_taylor(x)
local y = 1 - erfc(x) -- assume good erfc
—>y=.000000001128 on HP50G
return h==0 and y or y + 2/sqrt(pi) * exp(-x*x) * h
h=10^(-9)-8.86019642276E-10
But X small, x² very small, exp x² = 1 & exp (-x²) =1
—> y=y+2/sqrt(pi)*0*0=y+0=y

[Albert Chan]

Hi, Gil

It is just taylor series: y(x+h) ≈ y(x) + y'(x) * h

function ierf_taylor(x)
y = ierf(x)
x = erf(y)
dx = erf(y)' dy
dy/dx = 1/erf(y)' = sqrt(pi)/2 * exp(y*y)
approx(y) = ierfc(1-x) = erf(1-(1-x)) = erf(x-1+1)
h = exact - approx = x - (x-1+1)

FYI, if x is small, h is usually not zero.

>x = pi/1e6
>x - (x-1+1)
-3.4641E-13

Since small x correspond to small y, dy/dx is very flat ≈ sqrt(pi)/2

ierf(ε) ≈ sqrt(pi)/2 * (ε + pi/12*ε^3 + 7/480*pi^2*ε^5)

function erf_taylor(x)
y = erf(x)
dy = erf(x)' dx
dy/dx = erf(x)' = 2/sqrt(pi) * exp(-x*x)
approx(y) = 1 - erfc(x), we need to get what x 'belongs' to this y
h = exact - approx = x - ierf_taylor(y)

When x is small, dy/dx again is very flat ≈ 2/sqrt(pi)

erf(ε) ≈ 2/sqrt(pi) * (ε - ε^3/5 + ε^5/10)

[Albert Chan]

(11-14-2023 08:37 PM)Gil Wrote:  I try here
y = ierfc(1-x) -- assume good ierfc
With x=10^(-9) —> y=:x of Erfc: 8.86019642276E-10

This is very bad ierfc, only 3 digits accuracy, with exact input!
Perhaps UTPN is not so accurate after all?

ierf(1e-9) ≈ sqrt(pi)/2 * 1e-9 ≈ 8.86226925452758e-10 (all digits correct!)

Quote:function erf_taylor(x)
local y = 1 - erfc(x) -- assume good erfc
—>y=.000000001128 on HP50G
return h==0 and y or y + 2/sqrt(pi) * exp(-x*x) * h
h=10^(-9)-8.86019642276E-10
But X small, x² very small, exp x² = 1 & exp (-x²) =1
—> y=y+2/sqrt(pi)*0*0=y+0=y

Let's start from your y, but this time use a good ierf

lua> x, y = 1e-9, .000000001128
lua> dydx = 2/sqrt(pi) * exp(-x*x)
lua> h = x - ierf(y)
lua> y + dydx * h
1.1283791670955127e-09
lua> erf(x)
1.1283791670955127e-09

[Gil]

If x=10^(-9)
I try to improve the value found with my HP50G :
ERF(10^(-9)) : 000000001128.

Now with that value, .000000001128,
I calculate the inverse iERF(000000001128) :9.9911511335E-10.

It means that ERF (9.9911511335E-10) gives also 000000001128.

With Taylor
ERF (9.9911511335E-10 +8.8488665E-13)
= .000000001128 +f'×h, with h=8.8488665E-13

Correct?

And f'= 2/sqrt(pi)*exp(-x²).
But as mentioned in on HP50G, exp (very small) =1
And exp (-very small) = 1/1 = Constant on HP50G.
So that f'=2/sqrt(pi)*exp(-x²)
=2/sqrt(pi)*1 = 1.12837916709.
And f'×h = 1.12837916709 × 8.8488665E-13
= 9.98487661096E-13

And new found ERF =
000000001128 +f'×h=
000000001128 + 9.98487661096E-13=
1.12899848766E-9... ≠ 1.128379167!

I suppose that I understood all wrong.

Could you help me point the errors in order to get your very good solution value?

By the way, is there a way (if possible an easy one!), with the HP50G, to get an approximate value of
Exp(-1EE-18), which is less than 1?

Thanks for your comprehension and patience.

As you may have noticed, I am a dummy in maths and in many other topics.

Regards

[Albert Chan]

It seems a general solver to get from erfc to its inverse is a bad idea.

erfc(-∞) = 2
erfc(-5) ≈ 2
erfc(0) = 1
erfc(+5) ≈ 0
erfc(+∞) = 0

The curve is very flat on both end, it is hard to get accurate inverse, without function's derivatives.

f = a - erfc(x) = (a-1) + erf(x)
f' = 2/sqrt(pi) * exp(-x*x)
f''/f' = -2x

Halley: x = x - f/(f' - (f''/f')*(f/2)) = x - f/(f' + x*f)

Code:

function ierfc_halley(a, verbal)
    local x = ierf_rough(1-a) -- guess
    repeat
        local f = a - erfc(x)
        local h = -f/(2/sqrt(pi)*exp(-x*x)+x*f)
        if verbal then print(x, h) end
        x = x+h
    until x+h*h == x
    return x
end

function ierf_rough(x)  -- max rel err = 0.002
    local b, neg = 4.33, x<0
    local y = log(1-x*x)/2
    y = sqrt(sqrt(b*b+y*((2-pi)*b+y)) - (b+y))
    y = finite(y) or 6
    return neg and -y or y
end

Again, ierfc is only rely on accurate erfc, which we have.

lua> x = 1e-9
lua> y = ierfc_halley(1-x)
lua> y
8.862269003885489e-10
lua> ierfc(1-x)
8.862269003885489e-10

Perfect accuracy. Now continue from previously, to get ierf(x)

lua> dydx = sqrt(pi)/2 * exp(y*y)
lua> h = x - (x-1+1)
lua> y + dydx * h
8.86226925452758e-10
lua> ierf(x)
8.862269254527581e-10

[Gil]

1) Thanks for your fast reaction and your painstaking, Albert.

Will look at your explanations tomorrow.

2) To John : my apologises. Your test sign with IF condition was really necessary, for a zero value as input will indeed give an undermined sign ('?` on HP50G) : 2—nil for you! [/align]

3) Regarding the argument signs : I used induly that method for complex numbers, resulting in meaningless results —> the programs will have to be corrected.

[Gil]

I checked values given for UTPN (0, 1/2, x),
with x=8.81E-10, 8.82E-10,..., 8.99E-10.

All 12 significant digits always correct on the HP50G.

[Gil]

1) I corrected ERF and Erfc programs, in particular regarding complex numbers.

For ERFc(z), with z=0+bi, i.e. (0 b) in the program, I cheated and calculated
ERFc(0.0000000001 b), and replaced at the end the real part (about 0.9999...) by 1.

2) I checked the results with Wolfram Alpha and the HP PRIME.
- On the HP PRIME, approximate(erfc(5*i)), ie with no real part in the complex number input (argument), returns -82982738 79.76*i (a complex number with no real part).
- On my HP50G program ERFc(0 5) returns (1.,-8298273880.6 1 ), with 1 as the real part of that complex number.
- Whereas Wolfram Alpha returns 1 - 8.2982738806 7680 ×10^9 i.

Fazit
1) As no real part of the HP PRIME result is shown, that part has to be considered being equal to 0, and the PRIME result as wrong here (a bug?).
2) Imaginary part of the HP PRIME result is "worse" than my HP50G program (algorithm of the HP PRIME should be improved?).

3) Regarding Albert Chan observations
Still to be taken into account...

New corresponding codes

ERF
\\<< \"
1 Arg

.real 
x 
or
.complex 
(a b)
\" DROP DUP IM 0 ==
IF
THEN RE DUP 0 \\>=
IF
THEN 1
ELSE ABS -1
END
ELSE 1
END SWAP ERFc NEG 1 + *
\\>>

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*COS(2​*x*y)*SINH(k*y))))' \\->NUM * + - z 0 ==
IF
THEN C\\->R SWAP DROP 1 SWAP R\\->C
END
\\>>
END
\\>>

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

[Gil]

HP50: ERF(x), ERFc(x), iERFc(x) and iERFc(x), with x a real number
Now significand or mantissa with 11 to 12 accurate digits
(even for the inverse functions iERF and iERFc).

Note that using UTPN does not give always the full significant digits (or best solution).
The critical point of x lies then at about 0.114294.
Depending on the values of x, an integral was used accordingly (instead of UTPN).
—> For computational time, EMU48 is then highly recommended.

For ERF(z) and ERFc(z), with z a complex number (a b):
- the imaginary part of the resulting complex number should also give the correct 11-12 digits of the mantissa;
- the real part of the resulting complex number, however,
may not show all the mantissa digits (only some incomplete digits), depending on the case.
This issue still should be fixed up.

[Gil]

New version 6.

Better, but still not perfect.

1) Some input checking.

2) Improvements on final output and in the calculations.

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

4) In last version 5, I inadvertently overwrote the original inverse_ERF.

With my apologises.

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

[Gil]

New version 8

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

1) Same input checking as in versions 6/7
+ other input checkings,insuring that we are well in presence of not to large arguments (for real, max (IxI) OK if <= 9; for complex z=a+ib, max (IaI) and max(IbI) OK if both <= 9);

1b) The above max number (argument) 9 can be changed in the help program called max.arg.

1c) The help programs are labelled/written with small letters. See also 1g.

1d) The names of the 3 help programs are: complex (to deal with complex numbers), max.arg (above mentioned) and euler (to deal with complex numbers, to try and get the most exact result, when exact result should be ± (1 0) or ± (0 1)).

1e) The help programs are internal programs.
They are not to be used manually, or changed,
with the exception or the help program max.arg (see 1b).

1f) Instances of results with help program euler, in exact or approximate mode, with a dot after the 2 in the first two example:
'e^(2.*i*pi)' euler —> 1., and not (1.,2.54352307471E-12)
'e^(1.5*i*pi)' euler (0.,-1.), and not (-3.08985769397E-12,-1.)
but 'e^(2.2*i*pi)' euler —> (.809016994372,.587785252296).

1g) The name of all the user fonctions are written with (some) capital letters.
See also 1c.

2) Improvements on output "outlook".

3) Inverse functions, named now INV¦ and INV¦c
(and no more iERF and iERFc , as small i has another meaning —> see 4) below).

4) Besides the 3 functions related to the error function:
- iERFc
- ERFi
- w.FADDEEVA-KRAMP,
inclusion now of the following 3 new functions:
- PROBIT
- ERFcx
- F_DAWSON

For definitions, see Wikipedia/French and Wolfram Alpha, or 10) for the codes.

5) Inverse_ERF(p), labelled now INV¦
Initial program will convert p (p>=0 or p<0) into IpI.
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), called now INV¦c
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.9065 (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

6b) ERF and ERFc
The programs take into account the fact that the integral calculation gives more accurate result when IxI<= 0.01 then making use of the built-in function UTPN.

7) Initially, all programs allowed ± infinity as argument, supposing that the imaginary part of the argument did not exist. Now, it's not any more the case, but the "infinity" tests still appear, though they are not any more used in the present program execution. That part regarding "infinity" argument should be carefully removed in a cleaner version (or, again, treated separately, taking into account the fact that the argument may be a complex number (a+ib) with a=±infinity and b≠ infinity, or a≠infinity and b= ± infinity, or both a=±infinity and b=±infinity).

8) As in previous versions, all functions take one single argument.


9) Now all functions admit as argument, besides a real number, a complex number in the form:
- (a b)
- 'a+ib'
- 'r*e^(i*theta)'
- 'r*EXP(i*theta)',
even INV¦ and INV¦c functions (programs).

[Gil]

New version 9

MAIN FUNCTIONS:
all labelled with some capitals.

All of them accept real or complex format.

Automatic limitation of inputs:
to insure good results (up to 9-10 digits).

Many help programs/sub routines:
all labelled with small letters.

Arg/input = ± infinity, ±i × infinity:
ERF and ERFi functions OK for these extreme cases.

Reverse functions of ERF and ERFc only accept, as argument,
real x or z=a+ib with b=0.

All other functions accept, as argument,
real x or z=a+ib with a, b free real.

Extreme cases for arg/input = ± infinity, ±i × infinity
for functions iERFc ERFcx FADDEEVA DAWSON:
still to be handled.

[Gil]

Version 9b

Now inverse functions of the error function ERF and ERFc,
called respectively INV¦ & INV¦c,
admit as argument p a complex number p = a+ib, b≠0.

For the argument p of INV¦(p):
IpI < 1.

For the argument p of INV¦c(p):
I1-pI < 1.

On that purpose, a new subroutine called 'inv¦z' was used.
Better here to keep all the 50 terms of the numerator and of the denominator.
To save memory, however, you can use (approximate) real numbers, instead of the (exact) integers, in the 2 {} lists of the above mentioned subroutine 'inv¦z'.

Regarding PROBIT function
A slight change in the program, but the argument of PROBIT function was kept to remain a real number p in [0 1].

[Gil]

Version 9d

Much more compact subroutine 'inv¦z' for inverse calculation of ERF & ERFc,
when argument z is a complex with abs(z) < 1 :
now, to save memory place, the coefficients in the formulae development used in 'inv¦z' derive also from a formulae.

Strongly recommended to use the EMU48 Android application, if you want to get the results within a couple of seconds.