Functional inverse of $z=1+w+\cdots+w^{n-1}$ Migrated to MO.

I am interested in the functional inverse of
$$
z=1+w+\cdots+w^{n-1},\quad w\geq0,\ n>1.
$$
This function is strictly increasing on $w\geq0$ and thus admits an inverse.

My attempt:
By Lagrange's theorem we may write the inverse:
$$
w(z)=a+\sum _{k=1}^\infty g_k\frac{(z-f(a))^k}{k!},
$$
with $f(w)=1+w+\cdots+w^{n-1}$ and
$$
g_k=\lim_{w\to a}\partial_w^{k-1}\left({\frac {w-a}{f(w)-f(a)}}\right)^k.
$$
Expanding $f(w)$ around $w=1$ gives $f(w)=nF(1,1-n;2;1-w)$, which is a hypergeometric function. Choosing $a=1$ we write
$$
\frac{f(w)-f(1)}{w-1}=-n\frac{F(1,1-n,2,1-w)-1}{1-w}=\frac{n(n-1)}{2}F(1,2-n;3;1-w).
$$
It follows
$$
w(z)=1+\sum _{k=1}^\infty a_k\frac{2^k}{n^k(n-1)^k}\frac{(z-n)^k}{k!},
$$
with
$$
a_k=\lim_{w\to 1}\partial_w^{k-1}\left(F(1,2-n;3;1-w)\right)^{-k}.
$$
Evaluating the limit for $a_k$ is certainly non-trivial. One thought was to use the Faà di Bruno formula but I am a little unclear on the details of this calculation and am curious if $a_k$ can be written in a "nice" form that does not involve Bell polynomials. Could someone please fill me in on the details of computing the $a_k$'s? Also, for which values of $z$ does this series converge?
I was able to write a quick one line code in Mathematica to compute the $a_k$'s, which may lend itself to finding a pattern:
a[k_] := Limit[D[Hypergeometric2F1[1, 2 - n, c, 1 - w]^-k, {w, k - 1}], w -> 1] /. c -> 3

$a_1$, $a_2$, and $a_3$ seem to factor nicely while $a_4$ does not.

Edit:
Using Faà di Bruno's formula I was able to write down an explicit form for the $a_k$'s giving a final solution of
$$
\bbox[5px,border:2px solid #C0A000]{%
w(z)=1+\frac{2(z-n)}{n(n-1)}+\sum_{k=2}^\infty\sum_{\ell=1}^{k-1}(-k)^{(\ell)}B_{k-1,\ell}\left(\left\{(-1)^m\tfrac{m!(2-n)_m}{(3)_m}\right\}_{m=1}^{k-\ell}\right)\frac{\left(\frac{2(z-n)}{n(n-1)}\right)^k}{k!},%
}
$$
where $(s)^{(n)}=\Gamma(s+1)/\Gamma(s-n+1)$ is the falling factorial, $(s)_n=\Gamma(s+n)/\Gamma(s)$ is the Pochhammer symbol, and $B_{n,k}$ is the partial Bell polynomial. It is interesting to note that this result also works for the more general case $\{n\in\Bbb R:n>1\}$.
All that I am still curious about is the radius of convergence for this series which I believe is $|z-n|<n-1$.
Here is Mathematica code to compare the exact function $w(z)$ to approximation obtained by truncating its series expansion:
a[k_, n_] := 
 Sum[FactorialPower[-k, l] BellY[k - 1, l, 
    Table[(-1)^m (m! Pochhammer[2 - n, m])/Pochhammer[3, m], {m, 1, 
      k - l}]], {l, 1, k - 1}]
g[z_, n_, K_] := 
 1 + (2 (z - n))/(n (n - 1)) + 
  Sum[a[k, n] ((2 (z - n))/(n (n - 1)))^k/k!, {k, 2, K}]
gAprx[z_, n_] := 
 Quiet[N[Solve[(1 - w^n)/(1 - w) - z == 0, w, PositiveReals][[1, 1, 
     2]]]]
P[n_] := Manipulate[
  Plot[{gAprx[z, n], g[z, n, m]}, {z, 1, 2 n - 1}], {m, 2, 20, 1}]
P[5]

 A: Let $f:\mathbb C \to \mathbb C$ be the polynomial $\sum_{k=0}^{n-1} w^k$. We are interested in the radius of convergence of the power series of the inverse $w=f^{-1}(z)$ about $z=n$ corresponding to $w=1$. It is known that the radius of convergence of a power series represents the distance in the complex plane from the expansion point to the nearest singularity of the function expanded. Therefore, we have to locate the singularities of $f^{-1}$. The inverse of an entire function can have two types of singularities: $(1)$ algebraic singularities at the images of its saddle points $(2)$ transcendental singularities at its finite asymptotic values at infinity. There is no curve flowing to infinity along which a polynomial function tends to a finite value, hence the only singularities of our $f^{-1}$ will be algebraic singularities at the images of the saddle points of $f$. In formula, the singularities of $w=f^{-1}(z)$ are $z_k=f(w_k)$ where $w_k$ is one of the $n-2$ roots of the polynomial $f'(w)$. For the statement about the singularities of the inverse of an entire function, see Section 3, Appendix D in the book Lectures on the General Theory of Integral Functions by Georges Valiron.
A: Here is a solution based on the answer here which works as seen here when looking at the alternate forms section of the link. The derivation:
$$z=\sum_{k=0}^{n-1}w^k=\frac{w^n-1}{w-1}=a$$
so
$$w^n-1=a(w-1)=aw-a\implies w^n-aw+a-1=0$$
Remember the series definition of the Inverse of the Regularized Beta function where the series extends the domain of the function. It is also the quantile for a distribution using the Regularized Beta function:
$$\text I_z(a,b)=\frac{\int_0^z t^{a-1}(1-t)^{b-1}dt}{\text B(a,b)}=y\implies z=\text I^{-1}_y(a,b)$$
and using special cases of the Regularized Beta function:
$$\text I_{1-abw}(2,-n)=0\iff (c+1)a^n b^n w^n-abnw+n-1=0\implies w^n-\frac{abn}{a^n b^n (c+1)}w+\frac{n-1}{a^nb^n(c+1)}\mathop=^\text{set}w^n-aw+a-1=0$$
therefore:
$$c=n a^{-n}b^{1-n}-1:w^n-aw+\frac{n-1}{bn}$$
and
$$b=\frac{n-1}{(a-1)n}$$
Combining all of these identities, we put the problem in terms of the Regularized Beta Function:
$$w^n-aw+a-1=0\iff\text I_{1-\frac{a(n-1)}{(a-1)n}w}(2,-n)+\frac{(n-1)n^n (a-1)^{n-1}}{(a(n-1))^{n}}-1=0$$
Using the definition of the inverse and plugging in $w=0$:
$$a=\sum_{k=0}^{n-1}w^k=\frac{w^n-1}{w-1}\implies \boxed{w= -\frac{(a-1)n}{a(n-1)}\left(\text I^{-1}_{1-\left(\frac{a-1}{n-1}\right)^{n-1}\left(\frac na\right)^n}(2,-n)-1\right)}$$
where:
$$\boxed{\text I^{-1}_z(2,b)=\frac{2(b-1)z}{3b(b+1)}+\sqrt{\frac{2z}{b(b+1)}}+\frac{(b-1)(11b-2)\left(\frac z{b(b+1)}\right)^\frac32}{18\sqrt 2} +\frac{(b-1)(43b^2-5b-2)}{135b^2(b+1)^2}+… }$$
The function works for more general $a,w,n\in\Bbb R^+$,but here is the $a=5,n=7$ case:
$$w^7+…+1=5\iff w^7-5 w+ 4=0\implies w=.885128…= \frac{14}{15}-\frac{14}{15}\text I^{-1}_\frac{4246373}{56953125}(2,-7)$$
Using the $x$th partial sum and $y$=estimation, the above sum evaluates to:
$$(0,0.9333…)$$
$$(1,0.94217…)$$
$$(2,0.8865…)$$
$$(3,0.88482…)$$
$$(4,0.885197…)$$
Technically, $\text I^{-1}_z(a,b)$ is only defined for $a,b>0,0\le z\le1$, but the series extends the function’s domain. This method is more experimental, but it still works. The goal of this answer is the closed form, not the approximation. Please correct me and give me feedback!
Here is a graphical demo of the solution.
A: Extended Comment:
For the case $n=2$, the radius of convergence for the inverse series $w(z)$ is infinite, which can be deduced from the fact that $a_k=0$ for all $k\geq 2$.
For $n=3,4,\dots$, the approach given by @Gary provides a means for finding the radius of convergence of the inverse series either exactly or numerically. Let $w_k$ represent any one of the $n-2$ solutions of $f^\prime(w)=0$. Evaluating the derivative of $f(w)=(1-w^n)/(1-w)$ and equating with zero we find
$$
\begin{align}
f^\prime(w_k)=0
&\implies(n-1)w_k^n-nw_k^{n-1}+1=0,\quad w_k\neq 1\\
&\implies w_k^n=\frac{w_k}{n-(n-1)w_k},\quad w\neq 1\\
\end{align}
$$
Substituting $w_k^n$ into $f$ then gives
$$
z_k:=f(w_k)=\frac{n}{n-(n-1)w_k}.
$$
Now, after a bit of algebra
$$
|z_k-n|=(n-1)\biggl|\left(1+\tfrac{w_k}{n(1-w_k)}\right)^{-1}\biggr|;
$$
thus, the radius of convergence for any $n$ can be found by computing the minimum $|z_k-n|$, that is
$$
r(n)=\min_{k\in\{1,\dots,n-2\}}(n-1)\biggl|\left(1+\tfrac{w_k}{n(1-w_k)}\right)^{-1}\biggr|.
$$
This can be implemented in Mathematica as
fprime[w_, n_] := (1 - w^n)/(1 - w)^2 - ((n w^(-1 + n))/(1 - w));
wk[n_] := Solve[fprime[w, n] == 0, w][[All, 1, 2]];
ROC[n_] := 
  ToRadicals[
   Min[FullSimplify[(n - 1) Abs[(1 + wk[n]/(n (1 - wk[n])))^-1]]]];

For $n=3,4,5$ we obtain exact solutions, namely,
$$
\begin{align}
n=3: &|z-3|<\frac{9}{4}\\
n=4: &|z-4|<4\sqrt{\frac{2}{3}}\\
n=5: &|z-5|<\frac{5}{8 \sqrt{\frac{3}{140-5\ 2^{2/3} \sqrt[3]{\frac{5}{3 \sqrt{6}-2}}+5^{2/3} \sqrt[3]{2
   \left(3 \sqrt{6}-2\right)}}}}\\
\end{align}
$$
Mathematica is not able to find closed forms for larger values of $n$; however, we can compute the radius of convergence numerically with the modified code
fprime[w_, n_] := (1 - w^n)/(1 - w)^2 - ((n w^(-1 + n))/(1 - w));
wk[n_] := NSolve[fprime[w, n] == 0, w][[All, 1, 2]];
ROC[n_] := 
  Min[FullSimplify[(n - 1) Abs[(1 + wk[n]/(n (1 - wk[n])))^-1]]];

This code also works for noninteger values of $n$, e.g.
ROC[10.7]=9.75823

