Tunable Derivative-Matching Expansions

1. Introduction

Many approximations can be seen as “matching” procedures [1]: interpolations match function values, (generalized) Fourier series and moment-problem solutions match weighted integrals and several expansions are based on matching derivatives. Besides the Taylor series, these include a number of other approximations (Padé approximants, Neumann series of Bessel functions, Lambert series, etc...). Thus, once the derivatives of a function at some point are known, one can choose the approximation from a finite set of existing possibilities. With each of them reproducing the function well in the neighborhood of the expansion point, we can select the one that describes the function best when the argument becomes larger. For example, it is well known that for some applications the Padé rational function provides a better description than the Taylor series.

We adopt a similar strategy: we construct derivative-matching approximations with tunable parameters, meaning that for each (well-defined, non-singular) parameter combination the resulting approximation matches derivatives at the expansion point. One can then use the parameters to tune the far-away behavior. We propose six formulas, with five free parameters each. In this way six infinite sets of adjustable approximations are built.

Alternatively, one can also interpret our results as new formulas for special values of the Bell polynomials, since these are necessary for our purposes.

2. Method

The method of construction is not new and proceeds through a substitution into the standard power series $x\to g\left(x\right)$, as described in [2]. Reassessing our results there, we realized that some of them possess parametric freedom and this inspired us to investigate further. Adopting the notation

$${\mathcal{A}}^f\approx f\iff {\displaystyle \frac{d^n}{d{x}^n}}{\mathcal{A}}^f\left(x\right){{|}_{x=0}={\displaystyle \frac{d^n}{d{x}^n}}f\left(x\right)|}_{x=0},n\in {\mathbb{N}}_0$$

for an approximation ${\mathcal{A}}^f$, we assume, without loss of generality, that the expansion is performed at $x=0$. Writing ${\mathcal{A}}^f$ as

$${\mathcal{A}}^f\left(x\right)={\displaystyle \sum _{n=0}^{\infty }}{a}_n{\left[g\left(x\right)\right]}^n,g\left(0\right)=0,g^{\prime }\left(0\right)\ne 0,$$

(1)

straightforward manipulations and Faà di Bruno’s formula give

$$a_n={\displaystyle \frac{1}{n!}}{\displaystyle \frac{d^n}{d{x}^n}}f\left(g^{-1}\left(x\right)\right){|}_{x=0}={\displaystyle \frac{1}{n!}}{\displaystyle \sum _{k=0}^n}{d}_k^f{B}_{n,k}(d_1^{g^{-1}},d_2^{g^{-1}},\dots ,d_{n-k+1}^{g^{-1}}),$$

(2)

where $B_{n,k}$ are (incomplete) Bell polynomials and $d_n^h:={\displaystyle \frac{d^n}{d{x}^n}}h\left(x\right){|}_{x=0}$. The existence of a non-vanishing derivative of g at zero guarantees the existence of the inverse function $g^{-1}$ in some neighborhood of the same point.

Assuming we have a candidate function $g\left(x\right)$ to be used in (1), the simplest way to find $B_{n,k}:=B_{n,k}(d_1^{g^{-1}},d_2^{g^{-1}},\dots )$ is to study higher-order derivatives of ${\left(g^{-1}\right)}^k$. Indeed, the latter can be understood as $u\left(v\right(x\left)\right)$ where $v=g^{-1}$ and $u\left(x\right)=x^k$. Rewriting Faà di Bruno’s formula as

$$d_n^{u\circ v}={\displaystyle \sum _{j=0}^n}{d}_j^u{B}_{n,j}(d_1^v,d_2^v,\dots ),$$

one uses $d_j^u=k!{\delta }_{j,k}$ to get

$${\displaystyle \frac{d^n}{d{x}^n}}{\left(g^{-1}\left(x\right)\right)}^k{|}_{x=0}=k!{\displaystyle \sum _{j=0}^n}{\delta }_{j,k}{B}_{n,j}=k!B_{n,k}.$$

(3)

So, if one is able to derive an expression for $d_n^{{\left(g^{-1}\right)}^k}$ in a way that does not rely on the Bell polynomials, i.e., $d_n^{{\left(g^{-1}\right)}^k}=\mathrm{F}(n,k)$ then

$$B_{n,k}={\displaystyle \frac{1}{k!}}{d}_n^{{\left(g^{-1}\right)}^k}={\displaystyle \frac{1}{k!}}\mathrm{F}(n,k).$$

(4)

Thus, finding a formula for the nth derivative of ${\left(g^{-1}\right)}^k$ for some proposed g is one of the key concerns of this text. To implement parametric freedom, the free parameters are introduced into the definition of the function $g\left(x\right)=g(\left\{p_i\right\};x)$.

3. Properties

3.1. Existing Freedom

The Bell polynomials follow a transformation rule that reflects properties of differentiation

$$B_{n,k}\left(wq{c}_1,w{q}^2{c}_2,w{q}^3{c}_3,\dots \right)=q^n{w}^k{B}_{n,k}\left(c_1,c_2,c_3,\dots \right),$$

(5)

where q and w can be seen as parameters, which gives new possibilities for approximations. Indeed, if $f\left(x\right)\approx {\sum }_{n=0}^{\infty }{\displaystyle \frac{1}{n!}}\left({\sum }_{k=0}^n{d}_k^f{B}_{n,k}^{\circ }\right)g^n\left(x\right)$ then

$$f\left(x\right)\approx {\displaystyle \sum _{n=0}^{\infty }}{\displaystyle \frac{1}{n!}}\left({\displaystyle \sum _{k=0}^n}{d}_k^f{w}^k{B}_{n,k}^{\circ }\right)g^n\left({\displaystyle \frac{x}{w}}\right)$$

and

$$f\left(x\right)\approx {\displaystyle \sum _{n=0}^{\infty }}{\displaystyle \frac{1}{n!}}\left({\displaystyle \sum _{k=0}^n}{\displaystyle \frac{d_k^f}{q^n}}{B}_{n,k}^{\circ }\right){\left[qg\left(x\right)\right]}^n,$$

where $B_{n,k}^{\circ }:=B_{n,k}^{w=1,q=1}$. The latter expansion is trivial and does not represent a new approximation: the compensating scaling factor is introduced into the coefficients, but the summands remain identical. The former expansion is, however, non-trivial and we understand w as our first free parameter.

3.2. Convergence

Let I be the maximal open interval containing zero where g can be inverted. The question of convergence is then reduced to the convergence of the power series: the expansion $\sum a_n{g}^n\left(x\right)$ converges for all $x\in I$ such that $\left|g\left(x\right)\right|<R_y$, where $R_y$ is the radius of convergence of $\sum a_n{y}^n$. To determine $R_y$ one uses standard power-series methods. In general, the convergence to the function is lost outside I. If two different arguments $x_1\in I$ and $x_2\notin I$ give the same value $g\left(x_1\right)=g\left(x_2\right)$ then, obviously, the correct value is reconstructed for $x_1$ and the same value (which is generally incorrect) is obtained for $x_2$, since the series $\sum a_n{g}^n\left(x\right)$ cannot distinguish between $x_1$ and $x_2$.

Although convergence properties can be straightforwardly related to power series, the approximation properties are non-elementary: the expansion ${\mathcal{A}}^f=\sum a_n{g}^n$ of f may converge more rapidly and on a larger domain than the Taylor series of f. Consider a trivial example: with $g\left(x\right)=ln(1+x)$ and $f\left(x\right)=ln(1+x)$ the sum ${\mathcal{A}}^f$ approximates f exactly and very rapidly: just one term suffices. The Taylor series converge only at $(-1,1)$ and their convergence is slow (alternating sum). For $g\left(x\right)=ln(1+x)$, see the function $g_E$ from Definition 1 with $c=1$, $d=0$, $B=0$, $\alpha$=1 and $w=1$.

4. List

The following list of parametric functions $g_X$ consists of three pairs of inverse functions: any time a formula for $d_n^{{\left(g_X^{-1}\right)}^k}$ was established then also an expression for $d_n^{{\left(g_X\right)}^k}$ was found: the roles of $g_X$ and $g_X^{-1}$ can be swapped. It is unclear if this is just a coincidence or whether there is some deeper underlying reason. Parameters are included in the function names implicitly: $g_X\left(x\right)=g_X(p_1,p_2,p_3,p_4,p_5;x)$. Also we define $0^0:=1$.

Definition 1.

$$\begin{array}{cc}\hfill g_A\left(x\right)& ={\left({\displaystyle \frac{1+at+t^2}{1+bt+c{t}^2}}\right)}^{\alpha }-1,\mathrm{with}t={\displaystyle \frac{x}{w}}\hfill \\ \hfill g_B\left(x\right)& ={\displaystyle \frac{a-bt+\mathrm{sgn}(b-a)\sqrt{\left(b^2-4c\right)t^2+\left[4(c+1)-2ab\right]t+\left(a^2-4\right)}}{2\left[ct-1\right]}},\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill & \mathrm{with}t={\left(1+{\displaystyle \frac{x}{w}}\right)}^{1/\alpha }\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill g_C\left(x\right)& ={\displaystyle \frac{{\left[{(1+t)}^{\beta }+A\right]}^{\alpha }-{(1+A)}^{\alpha }}{1+B{\left[{(1+t)}^{\beta }+A\right]}^{\alpha }}},\mathrm{with}t={\displaystyle \frac{x}{w}}\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill g_D\left(x\right)& ={\left\{{\left[{\displaystyle \frac{t+{(1+A)}^{\alpha }}{1-Bt}}\right]}^{1/\alpha }-A\right\}}^{1/\beta }-1,\mathrm{with}t={\displaystyle \frac{x}{w}}\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill g_E\left(x\right)& ={\displaystyle \frac{{\left[ln(c+t)+d\right]}^{\alpha }-{\left[ln\left(c\right)+d\right]}^{\alpha }}{1+B{\left[ln(c+t)+d\right]}^{\alpha }}},\mathrm{with}t={\displaystyle \frac{x}{w}}\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill g_F\left(x\right)& =exp\left\{{\left[{\displaystyle \frac{t+{(lnc+d)}^{\alpha }}{1-Bt}}\right]}^{1/\alpha }-d\right\}-c,\mathrm{with}t={\displaystyle \frac{x}{w}}\hfill \end{array}$$

(11)

The Python implementation of generalized expansions using these functions is submitted as supplementary material with this text.

Proposition 2.

Assume $\alpha \ne 0$, $w\ne 0$ and $a\ne b$. Then $f\left(x\right)\approx {\sum }_{n=0}^{\infty }{a}_n{g}_A^n\left(x\right)$, where $a_n$ are defined by (2) with

$$B_{n,k}=w^k{\displaystyle \frac{n!}{k!}}{\displaystyle \sum _{m=k}^n}{\beta }_m{\displaystyle \sum _{l=0}^m}{(-1)}^{m-l}\left({\displaystyle \genfrac{}{}{0pt}{}{m}{l}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{l/\alpha }{n}}\right),$$

(12)

where

$$\begin{array}{cc}\hfill {\beta }_m& ={\displaystyle \frac{k}{m}}{\displaystyle \sum _{p=0}^{m-k}}{\displaystyle \sum _{j=0}^{\lfloor p/2\rfloor }}\left({\displaystyle \genfrac{}{}{0pt}{}{m}{p-j}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{p-j}{j}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{2m-k-p-1}{m-1}}\right)\hfill \\ \hfill & \times b^{p-2j}{c}^j{(a-b)}^{p+k-2m}{(c-1)}^{m-k-p}.\hfill \end{array}$$

Formula (6) has a significant limit case without the quadratic term in the numerator

$$g_A\left(x\right)=\underset{a,b,c,w\to \infty }{lim}{\left[{\displaystyle \frac{1+\frac{a}{w}x+\frac{1}{w^2}{x}^2}{1+\frac{b}{w}x+\frac{c}{w^2}{x}^2}}\right]}^{\alpha }-1,$$

where $a/w$, $b/w$ and $c/w^2$ tend to finite constants A, B and C. Thus

$$g_A\left(x\right)={\left[{\displaystyle \frac{1+Ax}{1+Bx+C{x}^2}}\right]}^{\alpha }-1={\left[{\displaystyle \frac{1+t}{1+K_{1}t+K_2{t}^2}}\right]}^{\alpha }-1,t={\displaystyle \frac{x}{w}},$$

(13)

with $w=1/A$, $K_1=B/A$, $K_2=C/A^2$, and one free parameter is lost. The corresponding limit of the formula for $C\ne 0$ affects only the ${\beta }_m$ term:

$$\begin{array}{cc}\hfill {\beta }_m^{\infty }& ={\displaystyle \frac{k}{m}}{\displaystyle \sum _{p=0}^{m-k}}{\displaystyle \sum _{j=0}^{\lfloor p/2\rfloor }}\left({\displaystyle \genfrac{}{}{0pt}{}{m}{p-j}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{p-j}{j}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{2m-k-p-1}{m-1}}\right)\hfill \\ \hfill & \times B^{p-2j}{C}^{j+m-k-p}{\left(A-B\right)}^{p+k-2m}.\hfill \end{array}$$

(14)

Taking $C=0$, one has

$${\beta }_m^{\infty ,C}=\left({\displaystyle \genfrac{}{}{0pt}{}{m-1}{k-1}}\right)B^{m-k}{\left(A-B\right)}^{-m}.$$

(15)

For the proofs, see the Appendix.

Proposition 3.

Assume $\alpha \ne 0$, $w\ne 0$, $a\ne b$ and $c\ne 1$. Then $f\left(x\right)\approx {\sum }_{n=0}^{\infty }{a}_n{g}_B^n\left(x\right)$, where $a_n$ are defined by (2) with

$$\begin{array}{cc}\hfill B_{n,k}& =w^k{\displaystyle \frac{n!}{k!}}{\displaystyle \sum _{j=1}^k}{(-1)}^{k-j}\left({\displaystyle \genfrac{}{}{0pt}{}{k}{j}}\right){\displaystyle \sum _{p=0}^n}{\displaystyle \sum _{q=0}^{\lfloor p/2\rfloor }}{a}^{p-2q}\left({\displaystyle \genfrac{}{}{0pt}{}{\alpha j}{p-q}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{p-q}{q}}\right)\hfill \\ \hfill & \times {\displaystyle \sum _{s=0}^{\lfloor (n-p)/2\rfloor }}{b}^{n-p-2s}{c}^s\left({\displaystyle \genfrac{}{}{0pt}{}{-\alpha j}{n-p-s}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{n-p-s}{s}}\right).\hfill \end{array}$$

(16)

For proof see Appendix.

Remark 4.

The case $c=1$ gives correct results if approached as the limit $c\to 1$.

Proposition 5.

Assume $\alpha \ne 0$, $\beta \ne 0$, $w\ne 0$, $A>-1$ and $B{(1+A)}^{\alpha }\ne -1$. Then $f\left(x\right)\approx {\sum }_{n=0}^{\infty }{a}_n{g}_C^n\left(x\right)$, where $a_n$ are defined by (2) with

$$\begin{array}{cc}\hfill B_{n,k}& =w^k{\displaystyle \frac{n!}{k!}}{\displaystyle \sum _{j=0}^k}\left({\displaystyle \genfrac{}{}{0pt}{}{k}{j}}\right){(-1)}^{k-j}{\displaystyle \sum _{p=0}^n}\left({\displaystyle \genfrac{}{}{0pt}{}{j/\beta }{p}}\right){\displaystyle \sum _{q=0}^p}\left({\displaystyle \genfrac{}{}{0pt}{}{p}{q}}\right){(-1)}^{p-q}\hfill \\ \hfill & \times {\displaystyle \sum _{r=0}^n}\left({\displaystyle \genfrac{}{}{0pt}{}{q/\alpha }{r}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{-q/\alpha }{n-r}}\right){(1+A)}^{p-\alpha r}{(-B)}^{n-r}.\hfill \end{array}$$

(17)

For the proof, see the Appendix.

Proposition 6.

Assume $\alpha \ne 0$, $\beta \ne 0$, $w\ne 0$, $A>-1$, $B{(1+A)}^{\alpha }\ne -1$ and $B\ne 0$. Then $f\left(x\right)\approx {\sum }_{n=0}^{\infty }{a}_n{g}_D^n\left(x\right)$, where $a_n$ are defined by (2) with

$$\begin{array}{cc}\hfill B_{n,k}& =w^k{\displaystyle \frac{n!}{k!}}{B}^{-k}{\displaystyle \sum _{j=0}^k}\left({\displaystyle \genfrac{}{}{0pt}{}{k}{j}}\right){(-1)}^j{\displaystyle \sum _{p=0}^n}\left({\displaystyle \genfrac{}{}{0pt}{}{-j}{p}}\right){\rho }^p{\displaystyle \sum _{q=0}^p}\left({\displaystyle \genfrac{}{}{0pt}{}{p}{q}}\right){(-1)}^{p-q}\hfill \\ \hfill & \times {\displaystyle \sum _{s=0}^n}\left({\displaystyle \genfrac{}{}{0pt}{}{\alpha q}{s}}\right)D^{-s}{\displaystyle \sum _{t=0}^s}\left({\displaystyle \genfrac{}{}{0pt}{}{s}{t}}\right){(-1)}^{s-t}\left({\displaystyle \genfrac{}{}{0pt}{}{\beta t}{n}}\right),\hfill \end{array}$$

(18)

where

$$D=1+A\mathrm{and}\rho ={\displaystyle \frac{B{D}^{\alpha }}{1+B{D}^{\alpha }}}.$$

For the proof, see the Appendix.

Remark 7.

The case $B=0$ gives correct results if approached as the limit $B\to 0$.

Proposition 8.

Assume $c>0$, $d+lnc>0$, $\alpha \ne 0$, $w\ne 0$ and $B{(d+lnc)}^{\alpha }\ne -1$. Then $f\left(x\right)\approx {\sum }_{n=0}^{\infty }{a}_n{g}_E^n\left(x\right)$, where $a_n$ are defined by (2) with

$$\begin{array}{cc}\hfill B_{n,k}& =w^k{\displaystyle \frac{n!}{k!}}{c}^k{\displaystyle \sum _{j=0}^k}\left({\displaystyle \genfrac{}{}{0pt}{}{k}{j}}\right){(-1)}^{k-j}{\displaystyle \sum _{p=0}^n}{\displaystyle \frac{j^p}{p!}}{\displaystyle \sum _{q=0}^p}\left({\displaystyle \genfrac{}{}{0pt}{}{p}{q}}\right){(-L_0)}^{p-q}\hfill \\ \hfill & \times {\displaystyle \sum _{r=0}^n}\left({\displaystyle \genfrac{}{}{0pt}{}{q/\alpha }{r}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{-q/\alpha }{n-r}}\right){(-B)}^{n-r}{L}_0^{q-\alpha r},\hfill \end{array}$$

(19)

where

$$L_0=lnc+d.$$

For the proof, see the Appendix.

Proposition 9.

Assume $c>0$, $d+lnc>0$, $\alpha \ne 0$, $w\ne 0$, $B{(d+lnc)}^{\alpha }\ne -1$ and $B\ne 0$. Then $f\left(x\right)\approx {\sum }_{n=0}^{\infty }{a}_n{g}_F^n\left(x\right)$, where $a_n$ are defined by (2) with

$$\begin{array}{cc}\hfill B_{n,k}& =w^k{\displaystyle \frac{B^{-k}}{k!c^n}}{\displaystyle \sum _{j=0}^k}\left({\displaystyle \genfrac{}{}{0pt}{}{k}{j}}\right){(-1)}^j{\displaystyle \sum _{p=0}^n}\left({\displaystyle \genfrac{}{}{0pt}{}{-j}{p}}\right){\left({\displaystyle \frac{B}{M}}\right)}^p{\displaystyle \sum _{q=0}^p}\left({\displaystyle \genfrac{}{}{0pt}{}{p}{q}}\right){(-L_0^{\alpha })}^{p-q}\hfill \\ \hfill & \times {\displaystyle \sum _{s=1}^n}s!\left({\displaystyle \genfrac{}{}{0pt}{}{\alpha q}{s}}\right)L_0^{\alpha q-s}s(n,s),\hfill \end{array}$$

(20)

where

$$L_0=lnc+d$$

and $s(n,s)$ are the Stirling numbers of the first kind.

For the proof, see the Appendix.

Remark 10.

The case $B=0$ gives correct results if approached as the limit $B\to 0$.

5. Case Studies

This chapter focuses on illustrative examples where proposed expansions are applied.

5.1. Increasing Convergence Rate and Domain of Power Series

In Section 3.2 we have demonstrated the convergence rate/domain enhancement using a somewhat degenerate (yet formally valid) case, in which the function being approximated was precisely the function used to build the series. Let us now provide a more “fair” example: we chose to approximate $f\left(x\right)=ln(1+x)$ using $g_C\left(x\right)$ with parameters $A=0$, $B=1$, $\alpha =1/2$, $\beta =1$ and $w=1$, thus obtaining

$$g_C\left(x\right)={\displaystyle \frac{\sqrt{1+x}-1}{\sqrt{1+x}+1}},$$

which is bijective on $[-1,+\infty )$. Besides elementary algebraic operations this formula uses only the square root which is evaluated very rapidly on a computer.

To investigate the behavior of the series

$$ln(1+x)\stackrel{?}{=}{\displaystyle \sum _{n=0}^{\infty }}{a}_n{g}_C^n\left(x\right)={\displaystyle \sum _{n=0}^{\infty }}{a}_n{\left({\displaystyle \frac{\sqrt{1+x}-1}{\sqrt{1+x}+1}}\right)}^n,$$

the inverse function is found to be

$$x={\left({\displaystyle \frac{1+g_C}{1-g_C}}\right)}^2-1.$$

It gives

$$\begin{array}{cc}\hfill ln(1+x)& =ln\left[{\left({\displaystyle \frac{1+g_C}{1-g_C}}\right)}^2\right]=2ln(1+g_C)-2ln(1-g_C).\hfill \end{array}$$

Using the standard expansion

$$ln(1+g_C)=g_C-{\displaystyle \frac{g_C^2}{2}}+{\displaystyle \frac{g_C^3}{3}}-{\displaystyle \frac{g_C^4}{4}}+\dots$$

and $|g_C\left(x\right)|<1$ for all $x>-1$, one sees that the expressions for $ln[1\pm g_C\left(x\right)]$ converge for all $x>-1$. One obtains

$$ln(1+x)=4\left(g_C+{\displaystyle \frac{g_C^3}{3}}+{\displaystyle \frac{g_C^5}{5}}+{\displaystyle \frac{g_C^7}{7}}\dots \right)={\displaystyle \sum _{n=0}^{\infty }}{a}_n{g}_C^n\left(x\right).$$

The summation no longer contains alternating signs, meaning that the expansion not only converges on the whole domain $(-1,+\infty )$ domain, but also converges more rapidly than the usual power expansion on $(-1,1)$. In addition, the convergence rate is enhanced by the absence of even powers.

In this studied example we used specific properties of the logarithm to find $a_{2n+1}=4/(2n+1)$, which is, of course, also reproduced by the formula (17). The formula is, however, much more general and covers many additional cases in which an increase in the convergence rate or domain may be achieved, and which are more difficult to analyze.

5.2. Reproducing Branch Points

Let us approximate the Lambert W function by constructing an expansion for its principal branch $W_0$ at $x=0$. We use $g_A$, formulas (13) and (15), with parameters $B=C=0$, $\alpha =1/2$ and $A=e$, where e is the base of the natural logarithm. Then

$$g_A\left(x\right)=\vartheta \sqrt{1+ex}-1,\vartheta =1.$$

(21)

This choice creates a branch point for $g_A$ situated at $x=-1/e$, i.e., at the exact same position as the branch point between $W_0$ and $W_{-1}$. Doing this, one hopes to construct an expansion that works beyond the branch cut and where the branch is chosen by selecting the appropriate branch of the square root in (21), i. e. $\vartheta =-1$ for $W_{-1}$.

Figure 1 shows the description of both $W_{0,-1}$ branches by $\sum a_n{g}_A^{\vartheta =1}$ and $\sum a_n{g}_A^{\vartheta =-1}$. One observes that the second branch is approximated as well which can be justified rigorously. Indeed, one can understand the two W branches as situated on the branch cut which starts at $x=-1/e$ and goes to $+\infty$. $W_0$ is then defined on its upper bank and $W_{-1}$ on the lower bank. Thus, $W_{-1}$ values are reached by transporting the $W_0$ values continuously around the branch point, not crossing the cut. It is known that the branch point is of the square-root type, so it can be removed by the substitution $t={(x+1/e)}^2-1/e$. Under this transformation, the upper bank remains unchanged (with a scaled argument), whereas the lower bank is “rotated” around $(-1/e,0)$ and identified with $(-\infty ,-1/e)$. Consequently, the branch point is removed. Specifically, let us define

$$R_W\left(x\right)=\left\{\begin{array}{cc}{W}_0\left[{(x+1/e)}^2-1/e\right]\hfill & \mathrm{if}x\ge -1/e\hfill \\ W_{-1}\left[{(x+1/e)}^2-1/e\right]\hfill & \mathrm{if}x<-1/e\hfill \end{array}\right..$$

From preceding arguments it follows that the function $R_W$ is analytic at $x=-1/e$ and the only (logarithmic) singularity of $R_W$ is situated at $x_L=-1/\sqrt{e}-1/e\approx -0.9744\dots$. By the transformation, it corresponds to the singularity of W on all Riemann sheets, $W_0$ excepted, at $x=0$. What is the effect of the argument transformation on $g_A$? One has

$$\begin{array}{cc}\hfill R_g\left(x\right)& =g_A\left(t\right)=\sqrt{e}\left(x+{\displaystyle \frac{1}{e}}\right)-1.\hfill \end{array}$$

We see that after the substitution the function used for the expansion becomes a simple linear function, $g_A\to R_g$, and the whole series $\sum a_n{g}_A^n$ becomes a derivative-matching power expansion $\sum a_n{R}_g^n$, which is the Taylor series. This converges to the nearest singularity which, after removing the $(-1/e,-1)$ branching, is the logarithmic singularity at $x_L$, see Figure 2. If $\sum a_n{R}_g^n\left(t\right)$ converges to the unified branches $W_{0,-1}\left(t\right)$, then, by back-substitution, $\sum a_n{g}_A^n{\left(x\right)|}_{\vartheta =\pm 1}$ converges to $W_0\left(x\right)$ and $W_{-1}\left(x\right)$, respectively. The back-substitution in fact gives a shifted Puiseux series.

Again, we have illustrated a feature of the approximations we propose on a simple example where the convergence can be easily demonstrated. Nevertheless, other parameter choices may also provide branch-point reconstruction (possibly with better convergence properties). For example, the numerical computations suggest that using $g_A$ with $b=c=0$, $\alpha =1/2$, $w=2$, $a=2e+1/\left(2e\right)$

$$g_A\left(x\right)=\vartheta \sqrt{1+\left(e+{\displaystyle \frac{1}{4e}}\right)x+{\displaystyle \frac{1}{4}}{x}^2}-1,\vartheta =1,$$

(22)

one gets series that converge to W in the neighborhood of the branch point for both branches (see Figure 3), although this may be more difficult to prove.

5.3. Matching Value

Matching the function value at a point $x_0\ne 0$, different from the expansion point, is possible only in special situations.

5.3.1. Scenario with $f\left(0\right)=0$ and $f\left(x_0\right)=0$

The generalized power series ${\sum }_{n=1}^N{a}_n{g}^n\left(x\right)$ takes the value zero at $x_0$ whenever $g\left(x_0\right)=0$, a condition that can be realized by suitable parameter tuning. However, such an $x_0$ necessarily lies outside the interval I containing 0, on which g is bijective (recall that $g\left(0\right)=0=g\left(x_0\right)$ while $g^{\prime }\left(0\right)\ne 0$). Consequently, in any neighborhood of $x_0$, the series will generally fail to converge to f. Nevertheless, this property remains useful, as it can still provide a rough approximation near specific points within the non-convergent region.

5.3.2. Scenario with $F\left(0\right)=F\left(X_0\right)=0$ and an Underlying Symmetry

In addition to $f\left(0\right)=f\left(x_0\right)=0$, let us also assume:

• The symmetry $g\left(x_1\right)=g\left(x_2\right)⟺f\left(x_1\right)=f\left(x_2\right)$;
• The convergence $\sum a_n{g}^n\left(x\right)\to f\left(x\right)$ for $x\in I$, where I is an interval containing zero on which g is bijective;
• The bound $f\left(x\right)\in J$ for all x, where J is the image of I, i.e., $J=f\left(I\right)$.

Proposition 11.

Under these strong conditions, we can conclude that

$$\sum _n{a}_n{g}^n\left(x\right)\to f\left(x\right)$$

for all x.

Proof. 

Indeed, for any $x_2$, we have $f\left(x_2\right)\in J=f\left(I\right)$, so there exists $x_1\in I$ such that $f\left(x_1\right)=f\left(x_2\right)$. Since $x_1\in I$, the assumed convergence gives $f\left(x_1\right)={\sum }_n{a}_n{g}^n\left(x_1\right).$ By the symmetry property, $f\left(x_1\right)=f\left(x_2\right)$ implies $g\left(x_1\right)=g\left(x_2\right)$. Consequently, ${\sum }_n{a}_n{g}^n\left(x_1\right)={\sum }_n{a}_n{g}^n\left(x_2\right)$, and hence

$$f\left(x_2\right)=f\left(x_1\right)=\sum _n{a}_n{g}^n\left(x_1\right)=\sum _n{a}_n{g}^n\left(x_2\right).$$

Since $x_2$ was arbitrary, the desired convergence holds for all x. □

Let us show how this property can be useful using an example that we investigate only numerically. We chose $f\left(x\right)=sin\left(\pi x\right)$ and $g_A$ with parameters $c=b=0$, $w=1$, $\alpha =1$ and $a=-1$, deliberately such that

$$g_A{\left(x\right)|}_{x=1}=\left(x^2-x\right){{|}_{x=1}=0=sin\left(\pi x\right)|}_{x=1},$$

i.e., the approximation reproduces the value of $sin\left(\pi x\right)$ at $x=1$. The symmetry of the functions with respect to $x=1/2$

$$g_A(1/2-x)=g_A(1/2+x),sin\left[\pi \left(1/2-x\right)\right]=sin\left[\pi \left(1/2+x\right)\right],$$

implies that if $\sum a_n{g}_A^n\left(x\right)\to sin\left(\pi x\right)$ for $0\le x\le 1/2$ then $\sum a_n{g}_A^n\left(x\right)\to sin\left(\pi x\right)$ also for $1/2\le x\le 1$. The imposed constraint ${\sum }_n{a}_n{g}_A^n\left(1\right)=0$ is expected to improve the approximation, for example when compared to the Taylor polynomial of the same degree. An objection can be made: because we add new information (about the symmetry), the improvement is not a surprise. This is only partially true. We may include the same information content also to an approximation based on a Taylor polynomial $T_N\left(x\right)={\sum }_{n=1}^N{\alpha }_n{x}^n$

$${\tau }_N\left(x\right)=\left\{\begin{array}{cc}{T}_N\left(x\right)\hfill & \mathrm{if}0\le x\le 1/2\hfill \\ T_N(1-x)\hfill & \mathrm{if}1/2<x\le 1\hfill \end{array}\right.,$$

(23)

thereby explicitly taking the symmetry into account. The results are shown in Figure 4. Clearly, the $g_A$-based approximation outperforms ${\tau }_N$ significantly. This can be understood: besides the symmetry property, our approximation is analytic at $x=1/2$ for any partial approximation ${\mathcal{A}}_N^f={\sum }_{n=1}^N{a}_n{g}_A^n$. We therefore believe that in a similar symmetry-backed situation the use of the value-matching procedure (value being zero) can be very useful and may considerably improve the convergence.

5.3.3. Optimizing Parameters

In previous paragraphs we have used specific parameters to demonstrate interesting features of new expansions without mentioning probably the most natural approach: tuning the parameters in an optimization procedure to best describe a function on a given interval. The aim would be to minimize the number of terms or the (computer) time to reach a desired precision. This ought to work: standard power expansions represent a single special case from the infinite set of parametric formulas we provide. There is no a priori reason why the global optimum should occur with the standard power series. Concerning computing time, they may win only if a small number of terms is taken into account, because our method requires computing $g\left(x\right)$ from x. For a larger number of terms, one can expect the benefits of using $g\left(x\right)$ (instead of x) to prevail. Let us also note that our expansions can be evaluated, like polynomials, by Horner’s method. Of course, for the generalized power series to be rapid, the coefficients need to be pre-computed, since they are computationally expensive. Thus a typical use case is to approximate a fixed function: compute the coefficients, store them, and use them repeatedly to evaluate the function for different arguments. We have made an effort to search for a mathematical function commonly implemented in numerical libraries (C, Java, Python, ...) for which our method could be competitive, at least on some reasonable interval. We did not succeed: the existing methods are fine-tuned over decades and use the power expansion (to be replaced by ours) very rarely; rather resorting to minimax polynomials and rational functions. Thus we believe our approach is more suited for specific, occasional situations where a solution to a problem is found as a non-standard (i.e., new, previously unknown) power series; it can then be enhanced.

Rather than providing a fully-tuned example let us focus on a specific optimization aspect: the parameter w is universal (present in all $g_X$) and computationally cheap. Tuning it may be a quick and routine way to improve the convergence. It does not appear deeply in formulas (12)-(20) and if these are used to compute $B_{n,k}^{\circ }:=B_{n,k}^{w=1}$ for $g\left(x\right)$ then $B_{n,k}=w^k{B}_{n,k}^{\circ }$, to be used with $g(x/w)$, see Section 3.1. It is therefore sufficient to remember a triangular array of numbers $B_{n,k}^{\circ }$ and search for an optimal w. The following example uses $g_F$ with $c=1$, $d=1$, $B=1/10$ and $\alpha =1$, i.e.,

$$g_F\left(x\right)=exp\left({\displaystyle \frac{11x}{10w-x}}\right)-1.$$

The series ${\sum }_{n=0}^N{a}_n{g}_F^n\left(x\right)$ with $a_n$ given by (20) is meant to approximate the Airy function Ai$\left(x\right)$ and we use the first six terms ($N=5$). In the optimization, we start with the reference value $w=1$ and tune it so that the approximation reproduces the first zero of Ai$\left(x\right)$, situated at $x=-2.33811\dots$. We find $w_{\mathrm{opt}.}=4.034\dots$ and present the results in Figure 5. We observe a very significant improvement and we see that the parameter tuning can hardly be compensated for by increasing the number of terms. We do not go beyond numerical observations in this example.

5.4. Manipulating Singularities

Assume f has a pole at $x_0$ and is analytic in its neighborhood $U=\{x:0<|x-x_0|<\epsilon \}$ for some $\epsilon >0$. It may then seem natural to approximate it using a function g that also has a pole at the same position. Unfortunately this does not work well.

Proposition 12.

Let us define ${\mathcal{A}}^f\left(x\right)={\sum }_{n=0}^{\infty }{a}_n{g}^n\left(x\right)$ and assume that the set $\{n:a_n\ne 0\}$ is infinite. Then there exists $x_1\in U$, $x_1\ne x_0$, such that ${\mathcal{A}}^f\left(x_1\right)\ne f\left(x_1\right)$.

Proof. 

By analyticity, we can write $f\left(x\right)={\sum }_{n=-N}^{\infty }{\alpha }_n{(x-x_0)}^n$ for $x\in U$ and separate the regular part $R_f\left(x\right)={\sum }_{n=0}^{\infty }{\alpha }_n{(x-x_0)}^n$ and the principal part $P_f\left(x\right)={\sum }_{n=1}^N{\alpha }_{-n}/{(x-x_0)}^n$. By definition $R_f$ is analytic (continuous) at $x=x_0$, so there exists ${\epsilon }_1\le \epsilon$ such that

$$|R_f\left(x\right)-{\alpha }_0|<{\displaystyle \frac{1}{2}}\mathrm{for}\mathrm{all}0<|x-x_0|<{\epsilon }_1.$$

This implies $|R_f\left(x\right)|<|{\alpha }_0|+{\displaystyle \frac{1}{2}}$. The principal part has a finite number of terms, and consequently its value is dominated by $1/{(x-x_0)}^N$ once x becomes close to $x_0$, i.e., there exists ${\epsilon }_2>0$ such that

$$|P_f\left(x\right)|>2|{\alpha }_0|+1\mathrm{for}\mathrm{all}0<|x-x_0|<{\epsilon }_2$$

Define ${\epsilon }_3=min\{{\epsilon }_1,{\epsilon }_2\}$ and $U_3=\{x:0<|x-x_0|<{\epsilon }_3\}$. Then, by the reverse triangle inequality

$$\left|f\left(x\right)\right|\ge |P_f\left(x\right)|-|R_f\left(x\right)|>(2|{\alpha }_0|+1)-\left(|{\alpha }_0|+{\displaystyle \frac{1}{2}}\right)=\left|{\alpha }_0\right|+{\displaystyle \frac{1}{2}}$$

for all $x\in U_3$. Since g has a pole at $x_0$ and infinitely many $a_n$ coefficients are nonzero, the function ${\mathcal{A}}^f$ has an essential singularity at $x_0$. If ${\mathcal{A}}^f$ does not converge for some $x\in U_3$, then the proof is complete. Otherwise, by Picard’s great theorem, ${\mathcal{A}}^f$ takes all complex values in $U_3$ with at most one exception $y_0$. Let $M=|{\alpha }_0|+1/2$ and define

$$S=\{x\in U_3:|{\mathcal{A}}^f\left(x\right)|<M,{\mathcal{A}}^f\left(x\right)\ne y_0\}.$$

Since ${\mathcal{A}}^f$ takes all values in $\mathbb{C}\setminus \left\{y_0\right\}$ infinitely many times in any neighborhood of $x_0$, the set S is nonempty. For any $x\in S$ we have

$$|{\mathcal{A}}^f\left(x\right)|<M=|{\alpha }_0|+{\displaystyle \frac{1}{2}}<\left|f\left(x\right)\right|$$

which implies ${\mathcal{A}}^f\left(x\right)\ne f\left(x\right)$. Thus ${\mathcal{A}}^f$ does not converge to f on S. □

If f has an essential singularity in $x_0$ then the convergence is possible. There are simple examples, see e.g., Section 5.1 of [2].

Although poles cannot be directly approximated by the poles of g, one may use poles of g to improve the convergence rate. For example, define g with a shape (concavity/convexity, inflection point) similar to f. Let us give an example with $f\left(x\right)=tan\left(x\right)$. We use

$$g_A\left(x\right)=i{\displaystyle \frac{x{x}_0}{x_0^2-x^2}},$$

(24)

which corresponds to parameters $\alpha =1$, $a=-1$, $b=0$, $c=1$ and $w=i{x}_0$ with $x_0=\pi (1+\sqrt{2})/2$. Here i is the imaginary unit and $\pm x_0$ are poles positions. One notices that $g_A\left(x\right)$ is purely imaginary on the real axis and its imaginary part is odd and similar to $tan\left(x\right)$, meaning it is increasing and antisymmetric with inflection point at the origin (Figure 6). By careful consideration of singularities of $tan\left(g^{-1}\right)$, the parameter $x_0$ is chosen so as to guarantee the convergence of the g power series on $(-\pi /2,\pi /2)$ (see Section 3.2). By combining the odd powers of g (even ones are missing) with imaginary coefficients, a real-valued function is constructed. Numerical inspection (Figure 7) suggests that the approximation converges more rapidly than the Taylor series of the same length. One notices that the formula (24) is a very simple rational function that is fast to evaluate on a computer.

6. Discussion, Summary, Outlook

The expansions we present are derived from the simple idea of a substitution into a power series; we do not see additional sources for our work that should imperatively be cited. After conducting extensive literature and internet research, we firmly believe that the presented results are original, nevertheless, we acknowledge the existence of results with some degree of similarity, which deserve to be mentioned. Many of them have a higher mathematical complexity and are often targeted at specific physics problems.

A well-known method for improving convergence of a power series is the Euler transform [3]

$$f\left(x\right)={\displaystyle \sum _{n=0}^{\infty }}{\alpha }_n{x}^n={\displaystyle \frac{1}{1-x}}{\displaystyle \sum _{n=0}^{\infty }}{a}_n{\left({\displaystyle \frac{x}{1-x}}\right)}^n,$$

(25)

where $a_n$ are computed as the nth order forward difference from coefficients ${\alpha }_j$. The method is well suited for summing alternating series and has generalization to negative indices [4] or to an approach dedicated to the summation of Rayleigh-Schrödinger series [5]. The formula (25) resembles what we present, but the factor $1/(1-x)$ reveals a difference.

Another mainstream approach for enhancing the convergence of series is Kummer’s transformation [6]. Aiming to speed up the convergence of $S={\sum }_{n=1}^{\infty }{\alpha }_n{x}^n$, one searches for a similar series $C={\sum }_{n=1}^{\infty }{\beta }_n{x}^n$ where C is known. Then S can be rewritten

$$S=\gamma C+{\displaystyle \sum _{n=1}^{\infty }}\left(1-\gamma {\displaystyle \frac{{\beta }_n}{{\alpha }_n}}\right){\alpha }_n{x}^n\equiv \gamma C+{\displaystyle \sum _{n=1}^{\infty }}{\lambda }_n{\alpha }_n{x}^n$$

with $\gamma ={lim}_{n\to \infty }{\alpha }_n/{\beta }_n$. If $\left\{{\beta }_n\right\}$ and $\left\{{\alpha }_n\right\}$ behave alike as n increases, then ${\lambda }_n$ tends to approach zero quickly.

A large number of publications stem from the seminal works of Aitken [7], Shanks [8] and Levin [9]. These authors introduced non-linear sequence transformations aimed at accelerating convergence by eliminating oscillations while preserving the underlying trend. Numerous other transformations have since been proposed [10,11,12], some of which incorporate free parameters. A nice overview is given in [13].

In fact, most of the methods mentioned above, as well as the method presented here, can be interpreted as conformal mappings, basically meaning substitution (or variable transformation), possibly with some terms re-arrangement. If done properly, it may move the point of evaluation from a position close to a singularity to a position far-away from the singularity, thus accelerating the convergence. Here, the essential texts to cite are by Van Dyke [14] and Guttmann [15].

An important approximation scheme is the Padé approximant which, by construction, is able to reproduce poles of the approximated function, thus possibly providing a good approximation even beyond this type of singularity. Here again, large generalizations exist: the authors of [16] propose a constructions where polynomials are replaced by general Chebyshev systems, thus possibly reproducing branch points and cuts. Another important application is the Padé-Borel summation, successfully implemented for example in quantum physics [17].

The latter citation is related to a topic many authors in physics and applied sciences work on: to give meaning to divergent series and to provide a mechanism for formally summing them up. Besides the Borel summation, one can mention order-dependent mapping methods [18], the variational perturbation theory [19] or more unconventional frameworks such as the self-similar approximation, a theory developed by Yukalov [20]. Nevertheless, these topics lie at the edge of the scope of this text.

We believe our construction surpasses previous approaches in several respects: the number and variety of both expansions and parameters, as well as their simplicity. The list of six functions we have presented is by no means exhaustive - one is free to invent his own g function, provided a formula for $d_n^{{\left(g^{-1}\right)}^k}$ can be derived. The same applies to the number and meaning of the parameters; we have focused on five-parameter formulas only for convenience, and one may attempt to construct a g with an even larger number of adjustable parameters.

Apart from regularizing divergent series, our approximations can be used, as demonstrated, for all the other purposes mentioned above. Yet the construction remains remarkably simple. In principle, Bell polynomials can be expressed by a formula, but this formula is formal rather than explicit. It contains a summation symbol where the sum runs over all sequences of indices subject to specific constraints; this does not qualify as an explicit formulation in our view. In contrast, the formulas we provide are truly explicit: although lengthy, they rely solely on single summations with explicit summation limits. Furthermore, the convergence of our expansions is easily determined via back-substitution into power series. None of the competing approaches combines all these features.

We cannot foresee all possible applications of our construction, yet, given the sheer size of the field dedicated to enhancing power expansions, we believe our method can make a meaningful contribution. It remains for each researcher or scientist to adapt the tools we propose to their own needs. Our experience shows that parameter tuning is not always easy; in preparing the examples for this paper, we have explored many suboptimal parameter regions and encountered problematic local minima during optimization. Still, we have always succeeded in finding an expansion that met our requirements. We hope this will likewise be the case in future applications.

As a future prospect, we plan to establish a diverse library of parametric functions similar to $g_X$ to provide a broader range of options.

Very generally, extensions of existing mathematical definitions or constructions can be useful; in the past, we [21,22] and others [23,24,25] have contributed to the concept of differentiation by integration, which generalizes the ordinary derivative (extends the set of differentiable functions, provides algorithms to differentiate noisy data or can be used for fractional derivative) and is parameterized by an infinite (but not arbitrary) set of kernel functions. We believe this is also a path to follow in the future: one may investigate useful extensions of existing mathematical concepts, which, in relation to the present text, may for example consist in considering expansions to other approximation forms (not only powers-based) and principles (not only derivative-matching).

Appendix A.1. Lagrange-Bürmann Formula

Let f be an analytic function and let h be its inverse near $x=0$. Then for any analytic function T

$$\left[x^n\right]T\left[h\left(x\right)\right]={\displaystyle \frac{1}{n}}\left[t^{n-1}\right]\left(T^{\prime }\left(t\right){\left[{\displaystyle \frac{t}{f\left(t\right)}}\right]}^n\right),$$

where $\left[x^n\right]F\left(x\right)$ denotes the coefficient of $x^n$ in the Taylor expansion of $F\left(x\right)$.

Appendix A.2. Proof of Formula ( uid11 )

Proof. 

We assume $\alpha \ne 0$ and $a\ne b$ (and $w=1$, $x/w\to x$). The inverse function to $g_A$, formula (6), is, in the neighborhood of $x=0$, the function $g_B$, formula (7), which can be verified by direct calculation. The latter is written as

$$g_B\left(x\right)={\displaystyle \frac{a-bu+\mathrm{sgn}(b-a)\sqrt{\Delta u}}{2\left[cu-1\right]}},$$

$$u={\left(1+x\right)}^{1/\alpha },\Delta u=\left(b^2-4c\right)u^2+\left[4(c+1)-2ab\right]u+\left(a^2-4\right).$$

One has

$$x=g_A\left[g_B\left(x\right)\right]={\left[{\displaystyle \frac{1+a{g}_B+g_B^2}{1+b{g}_B+c{g}_B^2}}\right]}^{\alpha }-1\Rightarrow {\left(x+1\right)}^{1/\alpha }={\displaystyle \frac{1+a{g}_B+g_B^2}{1+b{g}_B+c{g}_B^2}}.$$

Introducing

$$\omega :={\left(x+1\right)}^{1/\alpha }-1,\gamma :=a-b,\delta :=1-c$$

one has

$$\omega \left(g_B\right)={\displaystyle \frac{\gamma g_B+\delta g_B^2}{1+b{g}_B+c{g}_B^2}}.$$

Let us work in the variable $\omega$. We write $g_B^k\left(\omega \right)={\sum }_{m=k}^{\infty }{\beta }_m{\omega }^m$ where

$${\beta }_m={\displaystyle \frac{k}{m}}\left[t^{m-k}\right]{\left({\displaystyle \frac{1+bt+c{t}^2}{\gamma +\delta t}}\right)}^m.$$

(A1)

This follows from the Lagrange-Bürmann formula with $T\left(t\right)=t^k$, $h\left(\omega \right)=g_B\left(\omega \right)$:

$$\begin{array}{cc}\hfill {\beta }_m:=\left[{\omega }^m\right]{\left[g_B\left(\omega \right)\right]}^k& ={\displaystyle \frac{1}{m}}\left[t^{m-1}\right]\left(k{t}^{k-1}{\left[{\displaystyle \frac{t}{\omega {\left(g_B\right)}_{g_B=t}}}\right]}^m\right)\hfill \\ \hfill & ={\displaystyle \frac{k}{m}}\left[t^{m-k}\right]{\left({\displaystyle \frac{1+bt+c{t}^2}{\gamma +\delta t}}\right)}^m.\hfill \end{array}$$

We treat numerator and denominator in (A1) separately.

$$\begin{array}{cc}\hfill {\left(1+bt+c{t}^2\right)}^m& ={\left(1+X\right)}^m={\displaystyle \sum _{r=0}^m}\left({\displaystyle \genfrac{}{}{0pt}{}{m}{r}}\right)X^r={\displaystyle \sum _{r=0}^m}\left({\displaystyle \genfrac{}{}{0pt}{}{m}{r}}\right){\left(bt+c{t}^2\right)}^r\hfill \\ \hfill & ={\displaystyle \sum _{r=0}^m}\left({\displaystyle \genfrac{}{}{0pt}{}{m}{r}}\right){\displaystyle \sum _{j=0}^r}\left({\displaystyle \genfrac{}{}{0pt}{}{r}{j}}\right){\left(bt\right)}^{r-j}{\left(c{t}^2\right)}^j\hfill \\ \hfill & ={\displaystyle \sum _{r=0}^m}\left({\displaystyle \genfrac{}{}{0pt}{}{m}{r}}\right){\displaystyle \sum _{j=0}^r}\left({\displaystyle \genfrac{}{}{0pt}{}{r}{j}}\right)b^{r-j}{c}^j{t}^{r+j}{|}_{p=r+j}\hfill \\ \hfill & ={\displaystyle \sum _{p=0}^{2m}}{\displaystyle \sum _{j=0}^{⌊p/2⌋}}\left({\displaystyle \genfrac{}{}{0pt}{}{m}{p-j}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{p-j}{j}}\right)b^{p-2j}{c}^j{t}^p\hfill \end{array}$$

for any $m\in \mathbb{N}$. The original condition $r<m$ (i.e., $p-j<m$) is taken care of automatically by the binomial coefficient $\left({\displaystyle \genfrac{}{}{0pt}{}{m}{p-j}}\right)$. Next

$${(\gamma +\delta t)}^{-m}={\gamma }^{-m}{\displaystyle \sum _{q=0}^{\infty }}\left({\displaystyle \genfrac{}{}{0pt}{}{m+q-1}{q}}\right){\left(-{\displaystyle \frac{\delta }{\gamma }}\right)}^q{t}^q,$$

which follows from the generalized binomial theorem for negative exponents. One extracts the coefficient in front of the term $t^{m-k}$

$$\begin{array}{cc}\hfill \left[t^{m-k}\right]{\left({\displaystyle \frac{1+bt+c{t}^2}{\gamma +\delta t}}\right)}^m& =\left[t^{m-k}\right]{\displaystyle \sum _{p=0}^{2m}}{\displaystyle \sum _{j=0}^{⌊p/2⌋}}{\delta }_{p+q=m-k}\left({\displaystyle \genfrac{}{}{0pt}{}{m}{p-j}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{p-j}{j}}\right)\hfill \\ \hfill & \times \left({\displaystyle \genfrac{}{}{0pt}{}{m+q-1}{q}}\right)b^{p-2j}{c}^j{\gamma }^{-m}{\left(-{\displaystyle \frac{\delta }{\gamma }}\right)}^q{t}^{p+q}\hfill \\ \hfill & ={\displaystyle \sum _{p=0}^{2m}}{\displaystyle \sum _{j=0}^{⌊p/2⌋}}\left({\displaystyle \genfrac{}{}{0pt}{}{m}{p-j}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{p-j}{j}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{2m-k-p-1}{m-k-p}}\right)\hfill \\ \hfill & \times b^{p-2j}{c}^j{\gamma }^{p+k-2m}{(c-1)}^{m-k-p}.\hfill \end{array}$$

If $m-k-p<0$ (i.e., $p>m-k$) then the last binomial coefficient is zero so we can restrict the summation limit over p to $m-k$. Next, for the same binomial coefficient we use the identity $\left({\displaystyle \genfrac{}{}{0pt}{}{n}{r}}\right)=\left({\displaystyle \genfrac{}{}{0pt}{}{n}{n-r}}\right)$, thus getting

$$\begin{array}{cc}\hfill {\beta }_m& ={\displaystyle \frac{k}{m}}{\displaystyle \sum _{p=0}^{m-k}}{\displaystyle \sum _{j=0}^{⌊p/2⌋}}\left({\displaystyle \genfrac{}{}{0pt}{}{m}{p-j}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{p-j}{j}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{2m-k-p-1}{m-1}}\right)\hfill \\ \hfill & \times b^{p-2j}{c}^j{(a-b)}^{p+k-2m}{(c-1)}^{m-k-p}.\hfill \end{array}$$

(A2)

Let us recall our aim: we seek a formula for $d_n^{{\left(g_A^{-1}\right)}^k}$. We have $g_B^k\left(\omega \right)={\sum }_{m=k}^{\infty }{\beta }_m{\omega }^m$ and $\omega ={\left(x+1\right)}^{1/\alpha }-1$. Thus

$$\begin{array}{cc}\hfill d_n^{{\left(g_A^{-1}\right)}^k}& ={\displaystyle \frac{d^n}{d{x}^n}}{g}_B^k{\left(\omega \left(x\right)\right)|}_{x=0}\hfill \\ \hfill & ={\displaystyle \sum _{m=k}^{\infty }}{\beta }_m{\displaystyle \sum _{l=0}^m}\left({\displaystyle \genfrac{}{}{0pt}{}{m}{l}}\right){(-1)}^{m-l}{\displaystyle \frac{d^n}{d{x}^n}}{\left(x+1\right)}^{l/\alpha }{|}_{x=0}\hfill \\ \hfill & =n!{\displaystyle \sum _{m=k}^{\infty }}{\beta }_m{\displaystyle \sum _{l=0}^m}\left({\displaystyle \genfrac{}{}{0pt}{}{m}{l}}\right){(-1)}^{m-l}\left({\displaystyle \genfrac{}{}{0pt}{}{l/\alpha }{n}}\right)\hfill \end{array}$$

Using (4) we reproduce formula (12)

$$B_{n,k}={\displaystyle \frac{n!}{k!}}{\displaystyle \sum _{m=k}^{\infty }}{\beta }_m{\displaystyle \sum _{l=0}^m}{(-1)}^{m-l}\left({\displaystyle \genfrac{}{}{0pt}{}{m}{l}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{l/\alpha }{n}}\right),$$

with ${\beta }_m$ given by (A2). Special limit cases (14) and (15) are derived by applying the limit to the previous procedure. □

Appendix A.3. Proof of Formula ( uid16 )

Proof. 

We assume $\alpha \ne 0$, $a\ne b$ and $c\ne 1$ ($w=1$, $x/w\to x$). The inverse function to $g_B$, formula (7), is, in the neighborhood of $x=0$, the function $g_A$, formula (6), which can be verified by direct calculation. The latter is written as

$$g_A\left(x\right)={\left({\displaystyle \frac{1+ax+x^2}{1+bx+c{x}^2}}\right)}^{\alpha }-1.$$

So

$$\begin{array}{cc}\hfill {\left[g_A\left(x\right)\right]}^k& ={\left[{(1+ax+x^2)}^{\alpha }{(1+bx+c{x}^2)}^{-\alpha }-1\right]}^k\hfill \\ \hfill & ={\displaystyle \sum _{j=0}^k}\left({\displaystyle \genfrac{}{}{0pt}{}{k}{j}}\right){(-1)}^{k-j}{(1+ax+x^2)}^{\alpha j}{(1+bx+c{x}^2)}^{-\alpha j}\hfill \\ \hfill & \equiv {\displaystyle \sum _{j=1}^k}\left({\displaystyle \genfrac{}{}{0pt}{}{k}{j}}\right){(-1)}^{k-j}{F}_j\left(x\right),\hfill \end{array}$$

(A3)

where the lower limit can be set to $j=1$ because $j=0$ gives a constant which disappears in the differentiation. Now we search for $\left[x^n\right]F_j\left(x\right)=[d^n{F}_j\left(0\right)/d{x}^n]/n!$. Using the generalized binomial theorem we get

$$\begin{array}{cc}\hfill {(1+ax+x^2)}^{\alpha j}& ={\displaystyle \sum _{p=0}^{\infty }}\left({\displaystyle \genfrac{}{}{0pt}{}{\alpha j}{p}}\right){(ax+x^2)}^p\hfill \\ \hfill & ={\displaystyle \sum _{p=0}^{\infty }}\left({\displaystyle \genfrac{}{}{0pt}{}{\alpha j}{p}}\right){\displaystyle \sum _{q=0}^p}\left({\displaystyle \genfrac{}{}{0pt}{}{p}{q}}\right){\left(ax\right)}^{p-q}{\left(x^2\right)}^q\hfill \\ \hfill & ={\displaystyle \sum _{p=0}^{\infty }}{\displaystyle \sum _{q=0}^p}\left({\displaystyle \genfrac{}{}{0pt}{}{\alpha j}{p}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{p}{q}}\right)a^{p-q}{x}^{p+q}.\hfill \end{array}$$

To extract the coefficient in front of $x^r$ we introduce a new index $r=p+q$. Since $p>q$ one has $q\le ⌊r/2⌋$

$$\begin{array}{cc}\hfill \left[x^r\right]{(1+ax+x^2)}^{\alpha j}& =\left[x^r\right]{\displaystyle \sum _{q=0}^p}{\displaystyle \sum _{p=0}^{\infty }}{\delta }_{p+q=r}\left({\displaystyle \genfrac{}{}{0pt}{}{\alpha j}{p}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{p}{q}}\right)a^{p-q}{x}^{p+q}\hfill \\ \hfill & ={\displaystyle \sum _{q=0}^{⌊r/2⌋}}\left({\displaystyle \genfrac{}{}{0pt}{}{\alpha j}{r-q}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{r-q}{q}}\right)a^{r-2q}.\hfill \end{array}$$

Analogically

$$\left[x^t\right]{(1+bx+c{x}^2)}^{-\alpha j}={\displaystyle \sum _{s=0}^{⌊t/2⌋}}\left({\displaystyle \genfrac{}{}{0pt}{}{-\alpha j}{t-s}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{t-s}{s}}\right)b^{t-2s}{c}^s.$$

With $F_j\left(x\right)$ being a product of two terms, its $\left[x^n\right]$ factor can be expressed as a convolution

$$\begin{array}{cc}\hfill \left[x^n\right]F_j\left(x\right)& ={\displaystyle \sum _{p=0}^n}\left({\displaystyle \sum _{q=0}^{⌊p/2⌋}}\left({\displaystyle \genfrac{}{}{0pt}{}{\alpha j}{p-q}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{p-q}{q}}\right)a^{p-2q}\right)\hfill \\ \hfill & \times \left({\displaystyle \sum _{s=0}^{⌊(n-p)/2⌋}}\left({\displaystyle \genfrac{}{}{0pt}{}{-\alpha j}{n-p-s}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{n-p-s}{s}}\right)b^{n-p-2s}{c}^s\right).\hfill \end{array}$$

Substituting back to (A3) we get

$${\displaystyle \frac{d^n}{d{x}^n}}{\left[g_A\left(x\right)\right]}^k={\displaystyle \sum _{j=1}^k}\left({\displaystyle \genfrac{}{}{0pt}{}{k}{j}}\right){(-1)}^{k-j}{\displaystyle \frac{d^n}{d{x}^n}}{F}_j\left(x\right)=n!{\displaystyle \sum _{j=1}^k}\left({\displaystyle \genfrac{}{}{0pt}{}{k}{j}}\right){(-1)}^{k-j}\left[x^n\right]F_j\left(x\right).$$

This gives

$$\begin{array}{cc}\hfill B_{n,k}& ={\displaystyle \frac{n!}{k!}}{\displaystyle \sum _{j=1}^k}\left({\displaystyle \genfrac{}{}{0pt}{}{k}{j}}\right){(-1)}^{k-j}{\displaystyle \sum _{p=0}^n}\left({\displaystyle \sum _{q=0}^{⌊p/2⌋}}\left({\displaystyle \genfrac{}{}{0pt}{}{\alpha j}{p-q}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{p-q}{q}}\right)a^{p-2q}\right)\hfill \\ \hfill & \times \left({\displaystyle \sum _{s=0}^{⌊(n-p)/2⌋}}\left({\displaystyle \genfrac{}{}{0pt}{}{-\alpha j}{n-p-s}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{n-p-s}{s}}\right)b^{n-p-2s}{c}^s\right),\hfill \end{array}$$

which is identical to (16). □

Appendix A.4. Proof of Formula ( uid19 )

Proof. 

We assume $\alpha \ne 0$, $\beta \ne 0$, $A>-1$ and $B{(1+A)}^{\alpha }\ne -1$ ($w=1$, $x/w\to x$). The inverse function to $g_C$, formula (8), is, in the neighborhood of $x=0$, the function $g_D$, formula (9), which can be verified by direct calculation. The latter is written as

$$g_D\left(x\right)={\left\{{\left[{\displaystyle \frac{x+{(1+A)}^{\alpha }}{1-Bx}}\right]}^{1/\alpha }-A\right\}}^{1/\beta }-1.$$

We define $U\left(x\right)=g_D\left(x\right)+1$, so $g_D^k\left(x\right)={[U\left(x\right)-1]}^k$ and one has

$${\displaystyle \frac{d^n}{d{x}^n}}{\left[g_D\left(x\right)\right]}^k{{|}_{x=0}={\displaystyle \sum _{j=0}^k}\left({\displaystyle \genfrac{}{}{0pt}{}{k}{j}}\right){(-1)}^{k-j}{\displaystyle \frac{d^n}{d{x}^n}}{U}^j\left(x\right)|}_{x=0}.$$

Thus we search for $\left[x^n\right]U\left(x\right)$. We define

$$U\left(x\right)={\left[Z\left(x\right)-A\right]}^{1/\beta }\mathrm{with}Z\left(x\right):={\left[{\displaystyle \frac{x+{(1+A)}^{\alpha }}{1-Bx}}\right]}^{1/\alpha }.$$

We have $Z\left(0\right)=1+A$. So, introducing $D:=1+A$, we have

$$U^j\left(x\right)={\left[Z\left(x\right)-A\right]}^{j/\beta }={\left[\left(Z\left(x\right)-D\right)+\left(D-A\right)\right]}^{j/\beta }={\left[1+\left(Z\left(x\right)-D\right)\right]}^{j/\beta }.$$

The generalized binomial theorem gives

$$U^j\left(x\right)={\displaystyle \sum _{p=0}^{\infty }}\left({\displaystyle \genfrac{}{}{0pt}{}{j/\beta }{p}}\right){\left(Z\left(x\right)-D\right)}^p={\displaystyle \sum _{p=0}^{\infty }}\left({\displaystyle \genfrac{}{}{0pt}{}{j/\beta }{p}}\right){\displaystyle \sum _{q=0}^p}\left({\displaystyle \genfrac{}{}{0pt}{}{p}{q}}\right){(-D)}^{p-q}{Z}^q\left(x\right).$$

Next, let us focus on $Z^q\left(x\right)={(x+D^{\alpha })}^{q/\alpha }{(1-Bx)}^{-q/\alpha }$. We expand both parenthesis

$$\begin{array}{cc}\hfill {(x+D^{\alpha })}^{q/\alpha }& =D^q{\left(1+{\displaystyle \frac{x}{D^{\alpha }}}\right)}^{q/\alpha }={\displaystyle \sum _{r=0}^{\infty }}\left({\displaystyle \genfrac{}{}{0pt}{}{q/\alpha }{r}}\right)D^{q-\alpha r}{x}^r,\hfill \\ \hfill {(1-Bx)}^{-q/\alpha }& ={\displaystyle \sum _{s=0}^{\infty }}\left({\displaystyle \genfrac{}{}{0pt}{}{-q/\alpha }{s}}\right){(-B)}^s{x}^s.\hfill \end{array}$$

Again, we express $\left[x^n\right]Z^q\left(x\right)$ as convolution

$$\begin{array}{cc}\hfill \left[x^n\right]Z^q\left(x\right)& =\left[x^n\right]{\displaystyle \sum _{r=0}^{\infty }}{\displaystyle \sum _{s=0}^{\infty }}{\delta }_{r+s=n}\left({\displaystyle \genfrac{}{}{0pt}{}{q/\alpha }{r}}\right)D^{q-\alpha r}\left({\displaystyle \genfrac{}{}{0pt}{}{-q/\alpha }{s}}\right){(-B)}^s{x}^{r+s}\hfill \\ \hfill & ={\displaystyle \sum _{r=0}^n}\left({\displaystyle \genfrac{}{}{0pt}{}{q/\alpha }{r}}\right)D^{q-\alpha r}\left({\displaystyle \genfrac{}{}{0pt}{}{-q/\alpha }{n-r}}\right){(-B)}^{n-r}.\hfill \end{array}$$

Using ${\left[d^n{Z}^q\left(x\right)/d{x}^n\right]}_{x=0}=n!\left[x^n\right]Z^q\left(x\right)$ we get

$$\begin{array}{cc}\hfill {\displaystyle \frac{d^n}{d{x}^n}}{g}_D^k{\left(x\right)|}_{x=0}& ={\displaystyle \sum _{j=0}^k}\left({\displaystyle \genfrac{}{}{0pt}{}{k}{j}}\right){(-1)}^{k-j}{\displaystyle \frac{d^n}{d{x}^n}}{U}^j{\left(x\right)|}_{x=0}\hfill \\ \hfill & ={\displaystyle \sum _{j=0}^k}\left({\displaystyle \genfrac{}{}{0pt}{}{k}{j}}\right){(-1)}^{k-j}{\displaystyle \sum _{p=0}^{\infty }}\left({\displaystyle \genfrac{}{}{0pt}{}{j/\beta }{p}}\right){\displaystyle \sum _{q=0}^p}\left({\displaystyle \genfrac{}{}{0pt}{}{p}{q}}\right){(-D)}^{p-q}{\displaystyle \frac{d^n}{d{x}^n}}{Z}^q{\left(x\right)|}_{x=0}\hfill \\ \hfill & =n!{\displaystyle \sum _{j=0}^k}\left({\displaystyle \genfrac{}{}{0pt}{}{k}{j}}\right){(-1)}^{k-j}{\displaystyle \sum _{p=0}^{\infty }}\left({\displaystyle \genfrac{}{}{0pt}{}{j/\beta }{p}}\right){\displaystyle \sum _{q=0}^p}\left({\displaystyle \genfrac{}{}{0pt}{}{p}{q}}\right){(-1)}^{p-q}\hfill \\ \hfill & \times {\displaystyle \sum _{r=0}^n}\left({\displaystyle \genfrac{}{}{0pt}{}{q/\alpha }{r}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{-q/\alpha }{n-r}}\right)D^{p-\alpha r}{(-B)}^{n-r}.\hfill \end{array}$$

To complete the proof one needs to show that the upper summation index for p can be changed to n. If, for a fixed q, we define

$$F\left(q\right)={\displaystyle \sum _{r=0}^n}\left({\displaystyle \genfrac{}{}{0pt}{}{q/\alpha }{r}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{-q/\alpha }{n-r}}\right)D^{p-\alpha r}{(-B)}^{n-r}$$

then

$${\displaystyle \sum _{q=0}^p}\left({\displaystyle \genfrac{}{}{0pt}{}{p}{q}}\right){(-1)}^{p-q}F\left(q\right)={\Delta }^{p}F\left(0\right),$$

i.e., the sum over q represents pth finite forward difference. At the same time the expressions $\left({\displaystyle \genfrac{}{}{0pt}{}{q/\alpha }{r}}\right)$ and $\left({\displaystyle \genfrac{}{}{0pt}{}{-q/\alpha }{n-r}}\right)$ are, as function of q, polynomials of the order r and $n-r$ respectively. Thus their product has order n. Summing polynomials (index r) does not change the polynomial order. Therefore, if $p>n$, the expression becomes zero since any finite difference operator with a degree higher than is the degree of the polynomial on which it acts gives zero. So ${\sum }_{p=0}^{\infty }\to {\sum }_{p=0}^n$. Then, using $B_{n,k}=d_n^{{\left(g^{-1}\right)}^k}/k!$ we get formula (17). □

Appendix A.5. Proof of Formula ( uid21 )

Proof. 

We assume $\alpha \ne 0$, $\beta \ne 0$, $A>-1$, $B{(1+A)}^{\alpha }\ne -1$ and $B\ne 0$ ($w=1$, $x/w\to x$). The inverse function to $g_D$, formula (9), is, in the neighborhood of $x=0$, the function $g_C$, formula (8), which can be verified by direct calculation. The latter is written as

$$g_C\left(x\right)={\displaystyle \frac{{\left[{(1+x)}^{\beta }+A\right]}^{\alpha }-{(1+A)}^{\alpha }}{1+B{\left[{(1+x)}^{\beta }+A\right]}^{\alpha }}}.$$

We start by defining

$$W\left(x\right)={\left[{(1+x)}^{\beta }+A\right]}^{\alpha },D=1+A$$

giving

$$g_C\left(x\right)={\displaystyle \frac{W\left(x\right)-D^{\alpha }}{1+BW\left(x\right)}}={\displaystyle \frac{1}{B}}\left(1-{\displaystyle \frac{1+B{D}^{\alpha }}{1+BW\left(x\right)}}\right)\equiv {\displaystyle \frac{1}{B}}\left(1-{\displaystyle \frac{M}{1+BW\left(x\right)}}\right).$$

Then

$$g_C^k\left(x\right)=B^{-k}{\displaystyle \sum _{j=0}^k}\left({\displaystyle \genfrac{}{}{0pt}{}{k}{j}}\right){(-1)}^j{M}^j{\left[1+BW\left(x\right)\right]}^{-j}.$$

We remark that $W\left(0\right)=D^{\alpha }$ and we introduce $\Delta \left(x\right)$: $W\left(x\right)=D^{\alpha }+\Delta \left(x\right)$, where $\Delta \left(0\right)=0$. We get

$$1+BW\left(x\right)=1+B{D}^{\alpha }+B\Delta \left(x\right)=M+B\Delta \left(x\right)=M\left[1+{\displaystyle \frac{B}{M}}\Delta \left(x\right)\right].$$

Then the generalized binomial theorem is used

$${\left[1+{\displaystyle \frac{B}{M}}\Delta \left(x\right)\right]}^{-j}={\displaystyle \sum _{p=0}^{\infty }}\left({\displaystyle \genfrac{}{}{0pt}{}{-j}{p}}\right){\left({\displaystyle \frac{B}{M}}\right)}^p{\Delta }^p\left(x\right).$$

With $\Delta \left(0\right)=0$ one expands $\Delta \left(x\right)=0+\mathcal{O}\left(x\right)$. We are interested in the nth derivative at $x=0$, i.e., in the nth coefficient of the power expansion. The expression $\Delta \left(x\right)=\mathcal{O}\left(x\right)$ implies that ${\Delta }^p\left(x\right)$ contributes at order $\mathcal{O}\left(x^p\right)$, not earlier. We can safely truncate the series

$$\begin{array}{cc}\hfill {\left[1+BW\left(x\right)\right]}^{-j}& =M^{-j}{\displaystyle \sum _{p=0}^n}\left({\displaystyle \genfrac{}{}{0pt}{}{-j}{p}}\right){\left({\displaystyle \frac{B}{M}}\right)}^p{\Delta }^p\left(x\right)+\mathcal{O}\left(x^{n+1}\right)\hfill \\ \hfill & =M^{-j}{\displaystyle \sum _{p=0}^n}\left({\displaystyle \genfrac{}{}{0pt}{}{-j}{p}}\right){\left({\displaystyle \frac{B}{M}}\right)}^p{\displaystyle \sum _{q=0}^p}\left({\displaystyle \genfrac{}{}{0pt}{}{p}{q}}\right)W^q\left(x\right){\left(-D^{\alpha }\right)}^{p-q}+\mathcal{O}\left(x^{n+1}\right).\hfill \end{array}$$

Now we focus solely on the $W\left(x\right)$ term.

$$W\left(x\right)=D^{\alpha }{\left[1+{\displaystyle \frac{{(1+x)}^{\beta }-1}{D}}\right]}^{\alpha }\equiv D^{\alpha }{\left[1+T\left(x\right)\right]}^{\alpha },T\left(0\right)=0,$$

so

$$W^q\left(x\right)=D^{\alpha q}{\left[1+T\left(x\right)\right]}^{\alpha q}=D^{\alpha q}{\displaystyle \sum _{s=0}^{\infty }}\left({\displaystyle \genfrac{}{}{0pt}{}{\alpha q}{s}}\right)T^s\left(x\right).$$

Since $T\left(x\right)=\mathcal{O}\left(x\right)$, only terms $s=1,\dots ,n$ contribute to the nth coefficient. One has

$$\begin{array}{cc}\hfill \left[x^n\right]W^q\left(x\right)& =\left[x^n\right]D^{\alpha q}{\displaystyle \sum _{s=0}^n}\left({\displaystyle \genfrac{}{}{0pt}{}{\alpha q}{s}}\right)T^s\left(x\right)\hfill \\ \hfill & =D^{\alpha q}{\displaystyle \sum _{s=0}^n}\left({\displaystyle \genfrac{}{}{0pt}{}{\alpha q}{s}}\right)D^{-s}\left[x^n\right]{\left[{(1+x)}^{\beta }-1\right]}^s\hfill \\ \hfill & =D^{\alpha q}{\displaystyle \sum _{s=0}^n}\left({\displaystyle \genfrac{}{}{0pt}{}{\alpha q}{s}}\right)D^{-s}{\displaystyle \sum _{t=0}^s}\left({\displaystyle \genfrac{}{}{0pt}{}{s}{t}}\right){(-1)}^{s-t}\left[x^n\right]{(1+x)}^{\beta t}.\hfill \end{array}$$

With ${(1+x)}^{\beta t}={\sum }_{j=0}^{\infty }\left({\displaystyle \genfrac{}{}{0pt}{}{\beta t}{j}}\right)x^j$ obviously $\left[x^n\right]{(1+x)}^{\beta t}=\left({\displaystyle \genfrac{}{}{0pt}{}{\beta t}{n}}\right)$. This leads to

$$\left[x^n\right]W^q\left(x\right)=D^{\alpha q}{\displaystyle \sum _{s=0}^n}\left({\displaystyle \genfrac{}{}{0pt}{}{\alpha q}{s}}\right)D^{-s}{\displaystyle \sum _{t=0}^s}{(-1)}^{s-t}\left({\displaystyle \genfrac{}{}{0pt}{}{s}{t}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{\beta t}{n}}\right).$$

Now we put now everything together

$$\begin{array}{cc}\hfill {\displaystyle \frac{d^n}{d{x}^n}}{g}_C^k\left(x\right)& =\left[x^n\right]g_C^k\left(x\right)\hfill \\ \hfill & =B^{-k}{\displaystyle \sum _{j=0}^k}\left({\displaystyle \genfrac{}{}{0pt}{}{k}{j}}\right){(-1)}^j{M}^j\left[x^n\right]{\left[1+BW\left(x\right)\right]}^{-j}\hfill \\ \hfill & =B^{-k}{\displaystyle \sum _{j=0}^k}\left({\displaystyle \genfrac{}{}{0pt}{}{k}{j}}\right){(-1)}^j{M}^j{M}^{-j}{\displaystyle \sum _{p=0}^n}\left({\displaystyle \genfrac{}{}{0pt}{}{-j}{p}}\right){\left({\displaystyle \frac{B}{M}}\right)}^p\hfill \\ \hfill & \times {\displaystyle \sum _{q=0}^p}\left({\displaystyle \genfrac{}{}{0pt}{}{p}{q}}\right){\left(-D^{\alpha }\right)}^{p-q}\left[x^n\right]W^q\left(x\right).\hfill \end{array}$$

After back-substitutions, using $[d^n/d{x}^n]=n!\left[x^n\right]$ for $x=0$, formula (18) is obtained

$$\begin{array}{cc}\hfill B_{n,k}& ={\displaystyle \frac{n!}{k!}}{B}^{-k}{\displaystyle \sum _{j=0}^k}{(-1)}^j\left({\displaystyle \genfrac{}{}{0pt}{}{k}{j}}\right)\hfill \\ \hfill & \times {\displaystyle \sum _{p=0}^n}\left({\displaystyle \genfrac{}{}{0pt}{}{-j}{p}}\right){\left({\displaystyle \frac{B{D}^{\alpha }}{1+B{D}^{\alpha }}}\right)}^p{\displaystyle \sum _{q=0}^p}{(-1)}^{p-q}\left({\displaystyle \genfrac{}{}{0pt}{}{p}{q}}\right)\hfill \\ \hfill & \times {\displaystyle \sum _{s=0}^n}\left({\displaystyle \genfrac{}{}{0pt}{}{\alpha q}{s}}\right)D^{-s}{\displaystyle \sum _{t=0}^s}{(-1)}^{s-t}\left({\displaystyle \genfrac{}{}{0pt}{}{s}{t}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{\beta t}{n}}\right).\hfill \end{array}$$

□

Appendix A.6. Proof of Formula ( uid24 )

Proof. 

We assume $c>0$, $d+lnc>0$, $\alpha \ne 0$ and $B{(d+lnc)}^{\alpha }\ne -1$ ($w=1$, $x/w\to x$). The inverse function to $g_E$, formula (10), is, in the neighborhood of $x=0$, the function $g_F$, formula (11), which can be verified by direct calculation. The latter is written as

$$\begin{array}{cc}\hfill g_F\left(x\right)& =exp\left\{{\left[{\displaystyle \frac{x+{(lnc+d)}^{\alpha }}{1-Bx}}\right]}^{1/\alpha }-d\right\}-c=c\left[e^{Z\left(x\right)-L_0}-1\right]\hfill \end{array}$$

with

$$L_0=lnc+d\mathrm{and}Z\left(x\right)={\left({\displaystyle \frac{x+L_0^{\alpha }}{1-Bx}}\right)}^{1/\alpha }.$$

Then

$$g_F^k\left(x\right)=c^k{\left[e^{Z\left(x\right)-L_0}-1\right]}^k.$$

Setting $W\left(x\right)=Z\left(x\right)-L_0$, we get

$$\begin{array}{cc}\hfill g_F^k\left(x\right)& =c^k{\left[e^{W\left(x\right)}-1\right]}^k=c^k{\displaystyle \sum _{j=0}^k}\left({\displaystyle \genfrac{}{}{0pt}{}{k}{j}}\right){(-1)}^{k-j}{e}^{jW\left(x\right)}.\hfill \end{array}$$

We focus on the last term

$$e^{jW\left(x\right)}={\displaystyle \sum _{p=0}^{\infty }}{\displaystyle \frac{j^p}{p!}}{\left[Z\left(x\right)-L_0\right]}^p.$$

One has $Z\left(0\right)=L_0$ so $Z\left(x\right)-L_0=\mathcal{O}\left(x\right)$. Therefore ${\left[Z\left(x\right)-L_0\right]}^p$, when expanded, begins at $x^p$ or higher. Therefore

$$\begin{array}{cc}\hfill \left[x^n\right]e^{jW\left(x\right)}& =\left[x^n\right]{\displaystyle \sum _{p=0}^n}{\displaystyle \frac{j^p}{p!}}{\left[Z\left(x\right)-L_0\right]}^p\hfill \\ \hfill & ={\displaystyle \sum _{p=0}^n}{\displaystyle \frac{j^p}{p!}}{\displaystyle \sum _{q=0}^p}\left({\displaystyle \genfrac{}{}{0pt}{}{p}{q}}\right){(-L_0)}^{p-q}\left[x^n\right]Z^q\left(x\right).\hfill \end{array}$$

One has

$$\begin{array}{cc}\hfill Z^q\left(x\right)& ={\left({\displaystyle \frac{x+L_0^{\alpha }}{1-Bx}}\right)}^{q/\alpha }\hfill \\ \hfill & ={\left(x+L_0^{\alpha }\right)}^{q/\alpha }{\left(1-Bx\right)}^{-q/\alpha }\hfill \\ \hfill & ={\displaystyle \sum _{r=0}^{\infty }}\left({\displaystyle \genfrac{}{}{0pt}{}{q/\alpha }{r}}\right)x^r{L}_0^{q-\alpha r}{\displaystyle \sum _{s=0}^{\infty }}\left({\displaystyle \genfrac{}{}{0pt}{}{-q/\alpha }{s}}\right){\left(-B\right)}^s{x}^s\hfill \\ \hfill & ={\displaystyle \sum _{r=0}^{\infty }}{\displaystyle \sum _{s=0}^{\infty }}\left({\displaystyle \genfrac{}{}{0pt}{}{q/\alpha }{r}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{-q/\alpha }{s}}\right)L_0^{q-\alpha r}{\left(-B\right)}^s{x}^{r+s}.\hfill \end{array}$$

The coefficient in front of $x^n$ is given by the convolution $r+s=n$ ($s=n-r$)

$$\left[x^n\right]Z^q\left(x\right)={\displaystyle \sum _{r=0}^n}\left({\displaystyle \genfrac{}{}{0pt}{}{q/\alpha }{r}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{-q/\alpha }{n-r}}\right){\left(-B\right)}^{n-r}{L}_0^{q-\alpha r}.$$

This gives the final result

$$\begin{array}{cc}\hfill B_{n,k}& ={\displaystyle \frac{1}{k!}}{\displaystyle \frac{d^n}{d{x}^n}}{g}_F^k{\left(x\right)|}_{x=0}={\displaystyle \frac{n!}{k!}}\left[x^n\right]g_F^k\left(x\right)\hfill \\ \hfill & ={\displaystyle \frac{n!}{k!}}{c}^k{\displaystyle \sum _{j=0}^k}\left({\displaystyle \genfrac{}{}{0pt}{}{k}{j}}\right){(-1)}^{k-j}\left[x^n\right]{\displaystyle \sum _{p=0}^{\infty }}{\displaystyle \frac{j^p}{p!}}{\left[Z\left(x\right)-L_0\right]}^p\hfill \\ \hfill & ={\displaystyle \frac{n!}{k!}}{c}^k{\displaystyle \sum _{j=0}^k}\left({\displaystyle \genfrac{}{}{0pt}{}{k}{j}}\right){(-1)}^{k-j}{\displaystyle \sum _{p=0}^n}{\displaystyle \frac{j^p}{p!}}{\displaystyle \sum _{q=0}^p}\left({\displaystyle \genfrac{}{}{0pt}{}{p}{q}}\right){(-L_0)}^{p-q}\hfill \\ \hfill & \times {\displaystyle \sum _{r=0}^n}\left({\displaystyle \genfrac{}{}{0pt}{}{q/\alpha }{r}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{-q/\alpha }{n-r}}\right){\left(-B\right)}^{n-r}{L}_0^{q-\alpha r},\hfill \end{array}$$

which is identical to (19). □

Appendix A.7. Proof of Formula ( uid26 )

Proof. 

We assume $c>0$, $d+lnc>0$, $\alpha \ne 0$, $B{(d+lnc)}^{\alpha }\ne -1$ and $B\ne 0$ ($w=1$, $x/w\to x$). The inverse function to $g_F$, formula (11), is, in the neighborhood of $x=0$, the function $g_E$, formula (10), which can be verified by direct calculation. The latter is written as

$$g_E\left(x\right)={\displaystyle \frac{{\left[ln(c+x)+d\right]}^{\alpha }-{\left[ln\left(c\right)+d\right]}^{\alpha }}{1+B{\left[ln(c+x)+d\right]}^{\alpha }}}.$$

We define

$$V\left(x\right)={\left[ln(c+x)+d\right]}^{\alpha },L_0=lnc+d,M=1+B{L}_0^{\alpha }.$$

Then

$$\begin{array}{cc}\hfill g_E\left(x\right)& ={\displaystyle \frac{1}{B}}\left(1-{\displaystyle \frac{M}{1+BV\left(x\right)}}\right).\hfill \end{array}$$

So

$$g_E^k\left(x\right)=B^{-k}{\displaystyle \sum _{j=0}^k}\left({\displaystyle \genfrac{}{}{0pt}{}{k}{j}}\right){(-1)}^j{M}^j{\left[1+BV\left(x\right)\right]}^{-j}.$$

Writing $V\left(x\right)$ as $V\left(x\right)=L_0^{\alpha }+\Delta \left(x\right)$ with $\Delta \left(0\right)=0$ one has

$$\begin{array}{cc}\hfill 1+BV\left(x\right)& =M\left(1+{\displaystyle \frac{B}{M}}\Delta \left(x\right)\right).\hfill \end{array}$$

Raising to the power $-j$ gives

$${\left[1+BV\left(x\right)\right]}^{-j}=M^{-j}{\displaystyle \sum _{p=0}^{\infty }}\left({\displaystyle \genfrac{}{}{0pt}{}{-j}{p}}\right){\left({\displaystyle \frac{B}{M}}\right)}^p{\Delta }^p\left(x\right).$$

With $\Delta \left(x\right)=\mathcal{O}\left(x\right)$, the power expansion of ${\Delta }^p\left(x\right)$ starts at $x^p$ or higher. We are interested in extracting the nth coefficient, so we write

$$\begin{array}{cc}\hfill {\left[1+BV\left(x\right)\right]}^{-j}& =M^{-j}{\displaystyle \sum _{p=0}^n}\left({\displaystyle \genfrac{}{}{0pt}{}{-j}{p}}\right){\left({\displaystyle \frac{B}{M}}\right)}^p{\left[V\left(x\right)-L_0^{\alpha }\right]}^p+\mathcal{O}\left(x^{n+1}\right),\hfill \\ \hfill {\left[V\left(x\right)-L_0^{\alpha }\right]}^p& ={\displaystyle \sum _{q=0}^p}\left({\displaystyle \genfrac{}{}{0pt}{}{p}{q}}\right)V^q\left(x\right){(-L_0^{\alpha })}^{p-q}.\hfill \end{array}$$

This leads to

$$\left[x^n\right]{\left[1+BV\left(x\right)\right]}^{-j}=M^{-j}{\displaystyle \sum _{p=0}^n}\left({\displaystyle \genfrac{}{}{0pt}{}{-j}{p}}\right){\left({\displaystyle \frac{B}{M}}\right)}^p{\displaystyle \sum _{q=0}^p}\left({\displaystyle \genfrac{}{}{0pt}{}{p}{q}}\right){(-L_0^{\alpha })}^{p-q}\left[x^n\right]V^q\left(x\right).$$

Again, the generalized binomial theorem is used for the last term

$$\begin{array}{cc}\hfill \left[x^n\right]V^q\left(x\right)& =\left[x^n\right]{\left[ln(c+x)+d\right]}^{\alpha q}\hfill \\ \hfill & =\left[x^n\right]{\left[L_0+ln\left(1+{\displaystyle \frac{x}{c}}\right)\right]}^{\alpha q}\hfill \\ \hfill & ={\displaystyle \sum _{s=0}^{\infty }}\left({\displaystyle \genfrac{}{}{0pt}{}{\alpha q}{s}}\right)L_0^{\alpha q-s}\left[x^n\right]{\left[ln\left(1+{\displaystyle \frac{x}{c}}\right)\right]}^s.\hfill \end{array}$$

This can be related to the exponential generation function of the Stirling numbers of the first kind $s(n,k)$

$${\displaystyle \sum _{m=k}^{\infty }}s(m,k){\displaystyle \frac{x^m}{m!}}={\displaystyle \frac{{\left[ln\left(1+x\right)\right]}^k}{k!}}.$$

We have

$$\left[x^n\right]{\left[ln\left(1+{\displaystyle \frac{x}{c}}\right)\right]}^s=s!{\displaystyle \sum _{m=s}^{\infty }}s(m,s){\displaystyle \frac{1}{m!}}\left[x^n\right]{\displaystyle \frac{x^m}{c^m}}={\displaystyle \frac{s!}{n!c^n}}s(n,s).$$

Taking into the account the relation between power-expansion coefficients and derivatives, one reproduces the formula (20)

$$\begin{array}{cc}\hfill B_{n,k}& ={\displaystyle \frac{1}{k!}}{\displaystyle \frac{d^n}{d{x}^n}}{g}_E^k{\left(x\right)|}_{x=0}={\displaystyle \frac{n!}{k!}}\left[x^n\right]g_E^k\left(x\right)\hfill \\ \hfill & ={\displaystyle \frac{B^{-k}}{k!c^n}}{\displaystyle \sum _{j=0}^k}\left({\displaystyle \genfrac{}{}{0pt}{}{k}{j}}\right){(-1)}^j{\displaystyle \sum _{p=0}^n}\left({\displaystyle \genfrac{}{}{0pt}{}{-j}{p}}\right){\left({\displaystyle \frac{B}{M}}\right)}^p{\displaystyle \sum _{q=0}^p}\left({\displaystyle \genfrac{}{}{0pt}{}{p}{q}}\right){(-L_0^{\alpha })}^{p-q}\hfill \\ \hfill & \times {\displaystyle \sum _{s=1}^n}\left({\displaystyle \genfrac{}{}{0pt}{}{\alpha q}{s}}\right)s!L_0^{\alpha q-s}s(n,s).\hfill \end{array}$$

where the innermost sum was truncated thanks to relations $s(n,0)=0$ and $s(n,s)=0$ for $s>n$. □