The ancient pocket calculator HP32E computes $Q$, areas under the normal curve, with the "Algorithm 39" (same in FORTRAN) and the inverse $Q^{-1}$, quantile, using two different root finder. One is Newton's method, for the other I lack the know-why. I compared with Halley's method, with Householder's method, others, alas -- no match.
Details: By analysing the trace log of an emulated calculator (running the original firmware) I found, both functions, $Q$ and $Q^{-1}$, are divided in three sections using different formulas/procedures. Due to the symetry of the functions the outer sections use the same formula/procedure and adjust the result accordingly. For the middle section of $Q^{-1}$ Newton's method is used, while the algorithm for section I (and III) goes like this:
Setup of iteration
$m = \sqrt {-\mathrm{ln}\left (x^2\right )}$
$\displaystyle n = m - \frac{6,1}{\displaystyle 13 + m - \frac{94}{m + 8}}$,
what is (for $x \lt 0,1$) close to the asymptotic expansion $-\sqrt {\mathrm{ln} \frac{1}{x^2} - \mathrm{ln}\mathrm{ln} \frac{1}{x^2} - \mathrm{ln}(2\pi)}+o(1)$
$z_0 = n$, the initial guess, see note below after question
$h_0 = e^{\displaystyle n^2/2}$
The iteration
Do $i = 1$ until abs($k_i) <$ 1E-10
$r_i = r(z_{i-1})$ as in Q part 2, see also a. m. "Algorithm 39"
$s_i = 1/(z_{i-1} + b_2 + r_{i})$ with $b_2$ same as in Q part 2
$t_i = x \cdot {h_{i-1}}/{a_1}$ with $a_1 = 0,398942280444$ (what is almost $1/\sqrt{2\pi}$)
$k_i = z_{i-1}(t_i - s_i)$ what is $\displaystyle \sqrt{2\pi}\cdot h \cdot z \cdot (x-Q(z))$
$h_i = h_{i-1}/(1 + k_i)$
to apply the correction this way differs from Newton's method,
while for small $k$ it's quite close to the Taylor series.
$z_i = \sqrt{2 \mathrm{ln}(h_i)}$ what is $\sqrt{z^2-2 \mathrm{ln}(1+k)}$
loop/iterate/repeat
Question: What root finding method is this? Looking here and there for quite some time now, I found nothing similar yet.
Why I ask? I just like to grasp the math behind it.
Note: A detail nice to know, insignificant regarding my question. The HP32E sets the initial guess $z_0=\displaystyle n^2/2$ instead of $n$ as shown above. No doubt an error (IMO), which is ironed out by 1..2 more iterations.
Urged by the hint in the comment (links can rot away over time) here the formula "Q part 2" referenced above:

I documented it according the found procedure before I knew about its origin in "Algorithm 39".
Note, $t=x^2/2$
Some supplement: (mostly guesswork as I am not a mathematician)
My reasoning about the a. m. factor $\sqrt{2\pi}\cdot h \cdot z$ goes like this:
Comparing with Newton's method where the correction to $x_i$ is
$\frac{-\mathrm{f}(x_i)}{\mathrm{f'}(x_i)}$, the factor here must be deduced similarly.
In $\sqrt{z^2-2\,\mathrm{ln}(1+k)}$ the correction $k$ comprises $Q$ in some way, so this is the function which will with its derivative make the "funny factor" $\sqrt{2\pi}\cdot h \cdot z$
First, Q expressed by erf (error function):
$Q = \frac{\mathrm{erf}\left(\frac{x}{\sqrt{2}}\right)+1}{2}$
The iterated function, with (regarding the derivative) $k$ simplified to $k = Q$:
$\sqrt{z^2-2\,\mathrm{ln}(1+Q)}$
The reciprocal of its derivative (I got using Reduce CAS)
$\frac{1}{2}e^{\frac{x^2}{2}}\sqrt{-2\ln\left(\frac{\mathrm{erf}\left(\frac{x}{\sqrt{2}}\right)+3}{2}\right)+z^2}
\sqrt{2\pi}\left(-\mathrm{erf}\left(\frac{x}{\sqrt{2}}\right)-3\right)$
Simplified by let {sqrt(z^2-2*log(1 + Q)) => z, e^(x^2/2) => h}; results in
$\sqrt{2\pi}\cdot h \cdot z \cdot\frac{\left(-\mathrm{erf}\left(\frac{x}{\sqrt{2}}\right)-3\right)}{2}$
so here shows the "funny factor", alas, next step, final amelioration using the initial "Q by erf"
results in:
$$-\sqrt{2\pi}\cdot h \cdot z \cdot \left(Q+1\right)$$
So with few simplifications I get a formula with wrong sign and with $Q+1$ instead of $Q$. In addition, my "proof" is based on a step of the procedure I found analysing the HP32E, forming it back and forth to find the factor found in another step, kind of circular reasoning, I doubt, even if it looks nice.
Yes, numbers on the HP32E have a 10-digit mantissa for the user, internaly the exponent is saved in another register what enables computing with 13 digits. Results for the user are rounded. – m-stgt Jan 24 '24 at 03:03