InverseErf series expansion

In the Mathematica 11.1 documentation, they have given the following series expansion of the inverse error function. My question is: how does Mathematica compute it? Does it do an approximation of this series? If yes, how many terms do they include? Thanks.

Series[InverseErf[z], {z, 0, 10}]

(* (Sqrt[π] z)/2 + 1/24 π^(3/2) z^3 + 7/960 π^(5/2) z^5 + (127 π^(7/2) z^7)/80640 +
(4369 π^(9/2) z^9)/11612160 + O[z]^11 *)

For inverse functions like InverseErf[], the usual approach for determining their series coefficients is Lagrangian inversion. (That is, what InverseSeries[] does.)

In particular, Mathematica supports the (partial) Bell polynomials that appear in the series coefficients for the inverse series as BellY[]. (See also the discussion in Charalambides's book.) Using this, along with the general formula for the derivative of Erf[] (formula 7.1.19 in Abramowitz and Stegun), we have

SetAttributes[inverseErfCoefficient, Listable];
inverseErfCoefficient = 0; inverseErfCoefficient = Sqrt[π]/2;
inverseErfCoefficient[m_Integer?Positive] :=
(Sqrt[π]/2)^m BellY[Table[{(m + k - 2)!, (-2)^k/(2 k) Sqrt[π]/Gamma[1 - k/2]},
{k, 2, m}]]/(m! (m - 1)!)

Test:

inverseErfCoefficient[Range[0, 7]]
{0, Sqrt[π]/2, 0, π^(3/2)/24, 0, (7 π^(5/2))/960, 0, (127 π^(7/2))/80640}

Sum[inverseErfCoefficient[k] x^k, {k, 0, 20}] - (InverseErf[x] + O[x]^21)
O[x]^21

As it turns out, the series coefficients for the Taylor expansion of the inverse error function at $$0$$ satisfy a simple nonlinear recurrence. The formulae used in the following are adapted from Carlitz and Dominici:

ck = 0; ck = 1;
ck[n_Integer?Positive] := ck[n] = 2 Sum[Binomial[n - 1, k] ck[k - 1] ck[n - k], {k, n - 1}]

SetAttributes[inverseErfCoefficient, Listable];
inverseErfCoefficient[m_Integer?NonNegative] := (Sqrt[π]/2)^m ck[m]/m!

This version should give the same results.

• An example of why your answers are among my favorites. +1
– ciao
Nov 5 '17 at 6:32
• Glad to know you like it! I try my best when it comes to special functions... Nov 5 '17 at 7:25
• Thanks @J.M., I tested the code . It answers my question. However for larger values of x, say x=.9, I have to keep k=61 which is high. For x=.99 value of k is really high and it takes more time to compute so I was looking at the alternative . I came across another approximate formula on wikipedia, nand tested , it works better. it goes like this: anotherErfInv[x_] := Sqrt[Sqrt[(2/(.147 Pi ) + Log[1 - x^2]/2)^2 - Log[1 - x^2]/.147] - (2/(.147 Pi) + Log[1 - x^2]/2)]. This works better. Nov 6 '17 at 15:33
• @nutan, two things you need to remember: 1. a Taylor series is usually only accurate for values near its expansion point; and 2. no polynomial has ever had a vertical asymptote. Mull on those two statements, and then think about what you were doing again when you were evaluating at $x\approx 1$ Nov 6 '17 at 19:58
• This paper may be of interest: ams.org/journals/mcom/1976-30-136/S0025-5718-1976-0421040-7/… Though free access to the body of the paper is given at the AMS archive, the contents of the 62-page microfiche supplement (comprising 88 tables of rational approximants) are not present there. Nov 6 '17 at 22:17