How get a good numerical approximation for exponential quantile function?

Question

I want to calculate the "scaled" quantiles for an exponential distribution and I have a function for the inverse CDF, $iCDFExp(p,a,b)$ (e.g. $a=0, b=1$). The argument $p$ is in my application close to $1.0$, like $0.999\dots$, and in this situation, I get a numerical overflow.

Actually, the argument $p$ is also the result of a calculation, on a normal distribution:
I generate a uniform random variable $u$, convert it to a standard normal variable, and scale it up with a factor $s$ (e.g. $s=4$). This way I get an "upscaled" normal variable, and next, I convert it to uniform by using the CDF of the normal distribution. This way I get $p$, which I pass to the exponential inverse CDF:

$y=iCDFExp(CDFNorm(iCDFNorm(u)*s),a,b)$

Already for moderate $u$, like $0.99$, I get quite large normal variable, like $N=iCDFnorm(u)=2.5$, and with upscaling even $N*s=10$, so the CDFnorm is very close to $1.0$.

Is there an option to approximate $CDFNorm(iCDFNorm(u)*s)$ for large $u$, or even the full term $iCDFExp(CDFNorm(iCDFNorm(u)*s),a,b)$?

With $64$-bits the direct method starts to fail at $N=8.3$, but I need correct values also up to approximately $15$.

From simple thinking in terms of PDF, I would expect that y should roughly follow $s*s$.

Among the Wiki approximations I think that Winitzki fits best. So I wrote two functions.
Function CDFnorm_Winitzki(x: Double) : Double;
Const a=0.140012;
Var erf, xu, y : Double;
Begin
xu:=x/sqrt(2);
y:=(4/Pi+a * xu * xu)/(1+a * xu * xu);
erf:=signum(xu) * sqrt(1-exp(-xu * xu * y));
CDFnorm_Winitzki:=0.5+0.5 * erf
end;

Function iCDFnorm_Winitzki(x: Double) : Double;
Const a=0.140012;
Var erfinv, xu, y, lnx : Double;
Begin
xu:=2 * x-1;
c:=2/Pi/a;
lnx:=ln(1-xu * xu);
y:=sqrt(sqr(c + lnx/2) - lnx/a) - (c + lnx/2);
erfinv:=signum(xu) * sqrt(y);
iCDFnorm_Winitzki:=erfinv * sqrt(2); end;
But still I get the overflow!

Answer 1

Quoting the Wikipedia page

"The quantile function of a distribution is the inverse of the cumulative distribution function. The quantile function of the standard normal distribution is called the probit function, and can be expressed in terms of the inverse error function"

$$\Phi^{-1}(p) = \sqrt2\operatorname{erf}^{-1}(2p - 1)\qquad \text{for} \qquad p\in(0,1)$$

If you are not too requiring, may be you could use what I wrote here that is to say $$\text{erf}(x)\approx \sqrt{1-\exp\Big(-\frac 4 {\pi}\,\frac{1+\alpha\, x^2}{1+\beta \,x^2}\,x^2 \Big)}$$ where $$\alpha=\frac{10-\pi ^2}{5 (\pi -3) \pi } \qquad \text{and} \qquad \beta=\frac{120-60 \pi +7 \pi ^2}{15 (\pi -3) \pi }$$ which can easily be inversed.

This is similar to the approximation proposed by Sergei Winitzki (have a look here