{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "5a09a7f6",
   "metadata": {},
   "source": [
    "# Direct algorithm"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a002fd1d",
   "metadata": {},
   "source": [
    "We can calculate $[p/q](x)\\,$ using a direct method. Defining $N(x)$ and $D(x)$ as\n",
    "\\begin{equation*}\n",
    "N(x)=\\sum_{n=0}^{p} c_n\\,x^n\\,, \\qquad D(x)=\\sum_{n=0}^{q} b_n\\,x^n\\,,\n",
    "\\end{equation*} \n",
    "we have the following equivalence\n",
    "\\begin{equation*}\n",
    "f(x) - \\frac{\\sum_{n=0}^{p} c_n\\,x^n}{\\sum_{n=0}^{q} b_n\\,x^n}=\\mathcal{O}(x^{p+q+1}) \\quad \\Leftrightarrow \\quad\n",
    "\\left(\\sum_{n=0}^{+\\infty}a_n\\,x^n\\right) \\left(\\sum_{n=0}^{q} b_n\\,x^n\\right) - \\sum_{n=0}^{p} c_n\\,x^n = \\mathcal{O}(x^{p+q+1})\\,.\n",
    "\\end{equation*}   \n",
    "\n",
    "Matching the coefficients of $N(x)$ and $f(x)\\,D(x)\\,$ leads to a system of linear equations\n",
    "\n",
    "$$\\left\\{\\begin{array}{c c}\n",
    "\\begin{aligned}\n",
    "c_0&=a_0\\, b_0\\\\\n",
    "c_1&=a_1 \\,b_0 + a_0\\, b_1\\\\\n",
    "& \\ \\ \\ \\ \\ \\ \\ \\  \\ \\ \\ \\vdots \\ \\ \\ \\ \\ \\ \\ \\ \\vdots \\\\\n",
    "c_p&=a_p\\, b_0+ a_{p-1}\\,b_1+\\dots +a_{p-q}\\,b_q\\\\\n",
    "0&=a_{p+1}\\,b_0 + a_p\\,b_1+\\dots +a_{p-q+1}\\,b_q\\\\\n",
    "0&=a_{p+2}\\,b_0 + a_{p+1}\\,b_1+\\dots +a_{p-q+2}\\,b_q\\\\\n",
    "& \\ \\ \\ \\ \\ \\ \\ \\  \\ \\ \\ \\vdots \\ \\ \\ \\ \\ \\ \\ \\ \\vdots \\\\\n",
    "0&=a_{p+q}\\, b_0+ a_{p+q-1}\\,b_1+\\dots+ a_p\\,b_q \\,\\,\\\\\n",
    "\\end{aligned}\n",
    "\\end{array}\\right.$$\n",
    "\n",
    "with $a_n=0\\,,$  $\\forall n <0\\,.$ Selecting the last $q$ equations from the above system, we obtain the following matrix system  \t\t\n",
    "$$\\left(\\begin{array}{c c c c}\n",
    "a_{p-q+1}&a_{p-q+2}&\\dots \\dots&a_p\\\\\n",
    "a_{p-q+2}&a_{p-q+3}&\\dots \\dots&a_{p+1}\\\\\n",
    "\\vdots&\\vdots&\\ddots&\\vdots\\\\\n",
    "a_{p}&a_{p+1}&\\dots \\dots&a_{p+q-1}\\\\\n",
    "\\end{array}\\right)\n",
    "\\left(\\begin{array}{c}\n",
    "b_q\\\\\n",
    "b_{q-1}\\\\\n",
    "\\vdots\\\\\n",
    "b_1\\\\\n",
    "\\end{array}\\right)\n",
    "=\\left(\\begin{array}{c}\n",
    "-a_{p+1}\\\\\n",
    "-a_{p+2}\\\\\n",
    "\\vdots\\\\\n",
    "-a_{p+q}\n",
    "\\end{array}\\right)\\,,$$\n",
    "\n",
    "\n",
    "which can be expressed in a compact form\n",
    "\\begin{equation*}\\\n",
    "{\\bf A}\\,b=a \\,.\n",
    "\\end{equation*}\n",
    "When possible, there are several ways to solve the above system [1]. A classical approach is to use LU factorization, consisting of expressing ${\\bf A}={\\bf LU}$, it requires the order of $q^3$ operations to solve the system.\n",
    "Without loss of generality, assume that is possible to solve $${\\bf A}\\,b=a$$ and the denominator coefficients $b_0\\,,b_1\\,,b_2\\,, \\dots\\,,b_q\\,,$ are well defined. The remaining $p+1$ equations  \n",
    "\n",
    "$$\\begin{array}{c c}\n",
    "\\begin{aligned}\n",
    "c_0&=a_0\\,b_0\\,\\\\\n",
    "c_1&=a_1\\,b_0 + a_0\\, b_1\\\\\n",
    "c_2&=a_2\\,b_0 + a_1\\, b_1 + a_0\\,b_2\\\\\n",
    "\\ \\ \\  \\ \\ \\  \\ \\ & \\ \\ \\vdots\\\\\n",
    "c_p&=a_p\\, b_0+ a_{p-1}\\,b_1+\\dots +a_{p-q}\\,b_q\\\\\n",
    "\\end{aligned}\n",
    "\\end{array}$$\n",
    "gives the numerator coeficients $c_0\\,,c_1\\,,c_2\\,, \\dots\\,,c_p\\,$ values, and the direct algorithm ends. It requires the order of $q^3+pq\\,$ operations to derive an approximant.\n",
    "\n",
    "\n",
    "[1] Trefethen, L. N. and Bau, I. I. I. D. (1997). Numerical Linear Algebra. Society for Industrial and Applied\n",
    "Mathematics, Philadelphia."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "08d93b80",
   "metadata": {},
   "source": [
    "## Examples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "3d8d2faf",
   "metadata": {},
   "outputs": [],
   "source": [
    "import sympy as sp\n",
    "from padepy import direct_algorithm as da"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "ea180ae5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle e^{x}$"
      ],
      "text/plain": [
       "exp(x)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "var = sp.Symbol(\"x\")\n",
    "func = sp.exp(var)\n",
    "p = 3\n",
    "q = 1\n",
    "func"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "af3c0ea9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{\\frac{x^{3}}{24} + \\frac{x^{2}}{4} + \\frac{3 x}{4} + 1}{1 - \\frac{x}{4}}$"
      ],
      "text/plain": [
       "(x**3/24 + x**2/4 + 3*x/4 + 1)/(1 - x/4)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pade = da.pade(p, q, var, func)\n",
    "pade"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "4ee36a81",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{z + 1}{\\sqrt{z^{2} + 1}}$"
      ],
      "text/plain": [
       "(z + 1)/sqrt(z**2 + 1)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "var = sp.Symbol(\"z\")\n",
    "func = (var + 1) / sp.sqrt(var**2 + 1)\n",
    "p = 1\n",
    "q = 7\n",
    "func"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "142554f9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\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}$"
      ],
      "text/plain": [
       "(179*z/184 + 1)/(-5*z**7/128 + 147*z**6/1472 - 55*z**5/1472 - 31*z**4/368 - 15*z**3/368 + 97*z**2/184 - 5*z/184 + 1)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pade = da.pade(p, q, var, func)\n",
    "pade"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.9.5 64-bit",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.5"
  },
  "vscode": {
   "interpreter": {
    "hash": "62ffe33b4b70e70e2ad7eabb350bd852d57ace81263fde7960940ad59587f0df"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}