Baker’s algorithm

Henri Padé, in his Ph.D. thesis [1], organized the Padé approximants in a \(2\times2\) infinite table similar to next one.

\[\begin{split}\begin{array}{c| c c c c c c c c c c c c } [p / q] & 0 & 1 & 2 & 3 &\dots \\ \hline \\ 0 & [0/0] & [0/1] & [0/2] & [0/3] &\dots \\ \\ 1 & [1/0] & [1/1] & [1/2] & [1/3] &\dots \\ \\ 2 & [2/0] & [2/1] & [2/2] & [2/3] &\dots \\ \\ 3 & [3/0] & [3/1] & [3/2] & [3/3] &\dots \\ \vdots& \vdots &\vdots &\vdots&\vdots & \ddots \end{array}\end{split}\]

For different reasons we might be interested on particular blocks of the Padé table. On way of doing it is to use a recursive method. We now present Baker’s recursive algorithm [2].

Baker defined the sequences \begin{equation*} \begin{aligned} &\frac{\eta_{2j}(x)}{\theta_{2j}(x)}=[p-j/j\,]\,, \quad &\frac{\eta_{2j+1}(x)}{\theta_{2j+1}(x)}=[p-j-1/j\,]\, \end{aligned} \end{equation*} with \(j=0\,,1\,,\dots\,,p\,,\) and derived the recursive relations: \begin{equation*} \begin{aligned} &\frac{\eta_{2j}(x)}{\theta_{2j}(x)}=\frac{[\bar{\eta}_{2j-1}\,\eta_{2j-2}(x)-x\,\bar{\eta}_{2j-2}\,\eta_{2j-1}(x)]/\bar{\eta}_{2j-1}}{[\bar{\eta}_{2j-1}\,\theta_{2j-2}(x)-x\,\bar{\eta}_{2j-2}\,\theta_{2j-1}(x)]/\bar{\eta}_{2j-1}}\, \\ \\ &\frac{\eta_{2j+1}(x)}{\theta_{2j+1}(x)}=\frac{[\bar{\eta}_{2j}\,\eta_{2j-1}(x)-\bar{\eta}_{2j-1}\,\eta_{2j}(x)]/(\bar{\eta}_{2j}-\bar{\eta}_{2j-1})}{[\bar{\eta}_{2j}\,\theta_{2j-1}(x)-\bar{\eta}_{2j-1}\,\theta_{2j}(x)]/(\bar{\eta}_{2j}-\bar{\eta}_{2j-1})}\, \end{aligned} \end{equation*}

Where \(\bar{\eta}_{k}=\) is the coefficient of the higher power of \(\eta_{\,k}\) polynomial; otherwise, \(\bar{\eta}_{k}=0\). Both recursive relations are divided by a normalization constant \(\left(\theta_{k}(0)=1\,\right)\). It requires the order of \(p^2\) operations to derive an approximant. To start the recursive algorithm we input \begin{equation*} \begin{aligned} &\eta_{0}(x)=\sum_{n=0}^{p}a_n\, x^n\,, \quad &\theta_{0}(x)=1\,,\\ &\eta_{1}(x)=\sum_{n=0}^{p-1}a_n\, x^n \,, \quad &\theta_{1}(x)=1\,, \end{aligned} \end{equation*}

which follows the path in the Padé table illustrated below. The algorithm stops at \(j=p\,,\) with the calculation of \([0/p-1](x)\,\) and \([0/p](x)\,.\)

\[\begin{split}\begin{array}{| l |l|l|l|l|l|} \hline & & & & \ \ \ \ & \ \ \\ \hline & & & & &\\ \hline & & 5 \rightarrow& 6 \uparrow && \\ \hline & 3 \rightarrow& 4 \uparrow & & &\\ \hline 1 \rightarrow& 2 \uparrow & & & &\\ \hline 0 \uparrow & & & & &\\ \hline \end{array}\end{split}\]

For large values of \(p\) the method is clearly more efficient than the presented direct algorithm. Although, the direct algorithm might be better choice when the Padé table is seriously non-normal [2].

[1] Padé, H. (1892). Sur la représentation approchéee d’une fonction par des fractions rationnelles. Annales scientifiques de l’ École normale supérieure, 9(9):3-93.

[2] Baker, G. A. J. (1975). Essentials of Padé Approximants. Academic Press.

Examples

import sympy as sp
from padepy import baker_algorithm as ba

var = sp.Symbol("x")
func = sp.exp(var)
p = 2
q = 2
func

$\displaystyle e^{x}$

pade, path = ba.pade(p, q, var, func)
pade

Pade [2,2](x) stored at matrix index 4.

$\displaystyle \frac{\frac{x^{2}}{12} + \frac{x}{2} + 1}{\frac{x^{2}}{12} - \frac{x}{2} + 1}$

path

$\displaystyle \left[\begin{matrix}\frac{x^{4}}{24} + \frac{x^{3}}{6} + \frac{x^{2}}{2} + x + 1\\\frac{x^{3}}{6} + \frac{x^{2}}{2} + x + 1\\\frac{\frac{x^{3}}{24} + \frac{x^{2}}{4} + \frac{3 x}{4} + 1}{1 - \frac{x}{4}}\\\frac{\frac{x^{2}}{6} + \frac{2 x}{3} + 1}{1 - \frac{x}{3}}\\\frac{\frac{x^{2}}{12} + \frac{x}{2} + 1}{\frac{x^{2}}{12} - \frac{x}{2} + 1}\end{matrix}\right]$

var = sp.Symbol("z")
func = (var + 1) / sp.sqrt(var**2 + 1)
p = 3
q = 5
func

$\displaystyle \frac{z + 1}{\sqrt{z^{2} + 1}}$

pade, full_path = ba.pade(p, q, var, func, 0, False)
pade

Pade [3,5](z) stored at matrix index 10.

$\displaystyle \frac{\frac{147 z^{3}}{272} + \frac{5 z^{2}}{8} + \frac{147 z}{136} + 1}{- \frac{z^{5}}{64} + \frac{41 z^{4}}{272} + \frac{5 z^{3}}{136} + \frac{71 z^{2}}{68} + \frac{11 z}{136} + 1}$

full_path

$\displaystyle \left[\begin{matrix}\frac{35 z^{8}}{128} - \frac{5 z^{7}}{16} - \frac{5 z^{6}}{16} + \frac{3 z^{5}}{8} + \frac{3 z^{4}}{8} - \frac{z^{3}}{2} - \frac{z^{2}}{2} + z + 1\\- \frac{5 z^{7}}{16} - \frac{5 z^{6}}{16} + \frac{3 z^{5}}{8} + \frac{3 z^{4}}{8} - \frac{z^{3}}{2} - \frac{z^{2}}{2} + z + 1\\\frac{- \frac{75 z^{7}}{128} + \frac{z^{6}}{64} + \frac{45 z^{5}}{64} - \frac{z^{4}}{16} - \frac{15 z^{3}}{16} + \frac{3 z^{2}}{8} + \frac{15 z}{8} + 1}{\frac{7 z}{8} + 1}\\\frac{- \frac{11 z^{6}}{16} + \frac{7 z^{4}}{8} - \frac{3 z^{2}}{2} + 1}{1 - z}\\\frac{\frac{z^{6}}{64} - \frac{15 z^{5}}{352} - \frac{z^{4}}{16} + \frac{15 z^{3}}{44} + \frac{3 z^{2}}{8} + \frac{45 z}{44} + 1}{\frac{75 z^{2}}{88} + \frac{z}{44} + 1}\\\frac{- \frac{z^{5}}{24} - \frac{z^{4}}{24} + \frac{z^{3}}{3} + \frac{z^{2}}{3} + z + 1}{\frac{5 z^{2}}{6} + 1}\\\frac{- \frac{41 z^{5}}{704} + \frac{z^{4}}{16} + \frac{41 z^{3}}{88} + \frac{3 z^{2}}{4} + \frac{123 z}{88} + 1}{\frac{5 z^{3}}{16} + \frac{75 z^{2}}{88} + \frac{35 z}{88} + 1}\\\frac{- \frac{17 z^{4}}{56} - \frac{5 z^{2}}{7} + 1}{- \frac{11 z^{3}}{14} + \frac{11 z^{2}}{14} - z + 1}\\\frac{\frac{z^{4}}{16} + \frac{41 z^{3}}{68} + \frac{3 z^{2}}{4} + \frac{41 z}{34} + 1}{\frac{41 z^{4}}{272} + \frac{11 z^{3}}{68} + \frac{71 z^{2}}{68} + \frac{7 z}{34} + 1}\\\frac{\frac{z^{3}}{2} + \frac{z^{2}}{2} + z + 1}{\frac{z^{4}}{8} + z^{2} + 1}\\\frac{\frac{147 z^{3}}{272} + \frac{5 z^{2}}{8} + \frac{147 z}{136} + 1}{- \frac{z^{5}}{64} + \frac{41 z^{4}}{272} + \frac{5 z^{3}}{136} + \frac{71 z^{2}}{68} + \frac{11 z}{136} + 1}\\\frac{1 - \frac{23 z^{2}}{22}}{\frac{17 z^{5}}{88} - \frac{17 z^{4}}{88} - \frac{5 z^{3}}{11} + \frac{5 z^{2}}{11} - z + 1}\\\frac{\frac{5 z^{2}}{8} + \frac{147 z}{92} + 1}{\frac{147 z^{6}}{1472} - \frac{85 z^{5}}{736} - \frac{31 z^{4}}{368} + \frac{25 z^{3}}{92} + \frac{97 z^{2}}{184} + \frac{55 z}{92} + 1}\\\frac{z + 1}{\frac{z^{6}}{16} - \frac{z^{4}}{8} + \frac{z^{2}}{2} + 1}\\\frac{\frac{179 z}{184} + 1}{- \frac{5 z^{7}}{128} + \frac{147 z^{6}}{1472} - \frac{55 z^{5}}{1472} - \frac{31 z^{4}}{368} - \frac{15 z^{3}}{368} + \frac{97 z^{2}}{184} - \frac{5 z}{184} + 1}\\\frac{1}{- \frac{23 z^{7}}{16} + \frac{23 z^{6}}{16} - \frac{11 z^{5}}{8} + \frac{11 z^{4}}{8} - \frac{3 z^{3}}{2} + \frac{3 z^{2}}{2} - z + 1}\\\frac{1}{\frac{179 z^{8}}{128} - \frac{23 z^{7}}{16} + \frac{23 z^{6}}{16} - \frac{11 z^{5}}{8} + \frac{11 z^{4}}{8} - \frac{3 z^{3}}{2} + \frac{3 z^{2}}{2} - z + 1}\end{matrix}\right]$