Motivated by this now closed question, just for the art for art's sake, I wondered how efficiently we could approximate the inverse the function $$y=\frac{\big[\Gamma (1+x)\big]^a}{\Gamma (1+a x)}\qquad \text{with} \qquad 0 <y <1\qquad \text{and} \qquad a>1~~(a\in R)$$ Taking logarithms and using $$\log\big[\Gamma (1+\epsilon)\big]=-\gamma\, \epsilon +\sum_{n=2}^\infty (-1)^n\,\frac{\zeta(n)}{n}\, \epsilon^n$$ we have
$$-\log(y) =a\sum_{n=2}^\infty (-1)^{n} \,\frac{\zeta (n)}{n}\, (a^{n-1}-1)\, x^n$$ $$\color{red}{z}=\color{red}{-\frac{12 \log (y)}{\pi ^2 (a-1) a}}=x^2-\frac{4 (a+1) \zeta (3)}{\pi ^2}x^3+\frac{(a^2+a+1)\pi^2}{30} x^4-\frac{12 (a+1) \left(a^2+1\right) \zeta (5)}{5 \pi ^2}x^5+O\left(x^6\right)$$
Series reversion gives $$\color{blue}{x=z^{\frac 12}+\frac{2 (a+1) \zeta (3)}{\pi ^2}z+\left(\frac{10 (a+1)^2 \zeta (3)^2}{\pi ^4}-\frac{(a^2+a+1)\pi ^2}{60} \right)z^{\frac 32}+O(z^2)}$$
Some results for $a=\pi$ $$\left( \begin{array}{ccc} y & \text{approximation} & \text{solution} \\ 0.05 & 1.377403 & 1.359878 \\ 0.10 & 1.129186 & 1.118798 \\ 0.15 & 0.979501 & 0.972441 \\ 0.20 & 0.870317 & 0.865233 \\ 0.25 & 0.783314 & 0.779542 \\ 0.30 & 0.710272 & 0.707427 \\ 0.35 & 0.646767 & 0.644604 \\ 0.40 & 0.590124 & 0.588477 \\ 0.45 & 0.538587 & 0.537336 \\ 0.50 & 0.490918 & 0.489976 \\ 0.55 & 0.446195 & 0.445495 \\ 0.60 & 0.403684 & 0.403173 \\ 0.65 & 0.362757 & 0.362394 \\ 0.70 & 0.322831 & 0.322582 \\ 0.75 & 0.283299 & 0.283137 \\ 0.80 & 0.243435 & 0.243337 \\ 0.85 & 0.202211 & 0.202159 \\ 0.90 & 0.157825 & 0.157803 \\ 0.95 & 0.105843 & 0.105838 \end{array} \right)$$
which looks to be quite decent. In fact, it is better than the first iterate of Halley method with $x_0=\sqrt z$.
My question is
- Beside adding more terms to the expansion, will there be any other possible approach (independently of the accuracy) ?