InverseErf series expansion

Question

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 *)

Answer 1

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] = 0; inverseErfCoefficient[1] = 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] = 0; ck[1] = 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.