{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use of SymPy in SageMath."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1. Series Expansions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we use sympy in SageMath, then we may have to be explicit about the use of its functions.\n",
    "For example, the ``sqrt()`` function in SageMath is different from ``the sqrt()`` function in sympy.\n",
    "Consider the generating function for the Catalan numbers ``(1 - sqrt(4-x**2))/(2*x)``.\n",
    "With sympy we can define a generator for the terms in a Taylor series,\n",
    "but we must use the ``sqrt()`` function of sympy.  \n",
    "But it starts sooner, already with the type of variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<class 'sympy.core.symbol.Symbol'>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import sympy\n",
    "x = sympy.var('x')\n",
    "type(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Observe the difference with the var() that we are normally using in SageMath."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<class 'sage.symbolic.expression.Expression'>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y = var('y')\n",
    "type(y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Even if we use the sqrt() on a variable that is a sympy symbol, \n",
    "we get a symbolic expression of Sage on return."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<class 'sage.symbolic.expression.Expression'>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "type(sqrt(x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<class 'sage.symbolic.expression.Expression'>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "type(sqrt(y))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can make a generator object on the ``sqrt(1-x)``.  That works."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<generator object Expr.series.<locals>.<genexpr> at 0x7f587afddaf0>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tsqrtx = sympy.series(sqrt(1-x),x,0,n=None)\n",
    "tsqrtx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 1$"
      ],
      "text/plain": [
       "1"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "next(tsqrtx)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle - \\frac{x}{2}$"
      ],
      "text/plain": [
       "-x/2"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "next(tsqrtx)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[-x**2/8,\n",
       " -x**3/16,\n",
       " -5*x**4/128,\n",
       " -7*x**5/256,\n",
       " -21*x**6/1024,\n",
       " -33*x**7/2048,\n",
       " -429*x**8/32768,\n",
       " -715*x**9/65536]"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[next(tsqrtx) for _ in range(8)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But it will not work on ``sqrt(1-y)``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "tsqrty = sympy.series(sqrt(1-y),y,0,n=None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<class 'generator'>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "type(tsqrty)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\sqrt{1 - y}$"
      ],
      "text/plain": [
       "sqrt(1 - y)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "next(tsqrty)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Doing another ``next(tsqrty)`` triggers the ``StopIteration`` exception."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Although we can explitly use the sqrt() of sympy,\n",
    "it does make no difference, because ``y`` is a SageMath variable."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So for the generating function for the Catalan numbers, \n",
    "we must use a sympy symbol x and not a Sage variable y."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1,\n",
       " x,\n",
       " 2*x**2,\n",
       " 5*x**3,\n",
       " 14*x**4,\n",
       " 42*x**5,\n",
       " 132*x**6,\n",
       " 429*x**7,\n",
       " 1430*x**8,\n",
       " 4862*x**9,\n",
       " 16796*x**10]"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "g = (1 - sqrt(1-4*x))/(2*x)\n",
    "cg = sympy.series(g,x,0,n=None)\n",
    "[next(cg) for k in range(11)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "16796"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "catalan_number(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. Pattern Matching"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another feature that is a bit better in SymPy than in Sage is pattern matching,\n",
    "which can be useful when combined with substitution.\n",
    "To run this code we can also open a Terminal session, select Misc, and select python."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{y_: 10}"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "reset()\n",
    "from sympy import *\n",
    "x = Symbol('x')\n",
    "y = Wild('y')\n",
    "d = (10*x**3).match(y*x**3)\n",
    "d"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. Solving Recurrence Relations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "SymPy can solve recurrence relations.\n",
    "Consider the Fibonacci numbers,\n",
    "defined as $f(n) = f(n-1) + f(n-2)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle - f{\\left(n \\right)} - f{\\left(n + 1 \\right)} + f{\\left(n + 2 \\right)}$"
      ],
      "text/plain": [
       "-f(n) - f(n + 1) + f(n + 2)"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "reset()\n",
    "from sympy import Function, rsolve\n",
    "from sympy.abc import n\n",
    "f = Function('f')\n",
    "Fib = f(n+2) - f(n+1) - f(n)\n",
    "Fib"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To solve a recurrence relation, we call ``rsolve``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle C_{0} \\left(\\frac{1}{2} - \\frac{\\sqrt{5}}{2}\\right)^{n} + C_{1} \\left(\\frac{1}{2} + \\frac{\\sqrt{5}}{2}\\right)^{n}$"
      ],
      "text/plain": [
       "C0*(1/2 - sqrt(5)/2)**n + C1*(1/2 + sqrt(5)/2)**n"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "s = rsolve(Fib, f(n))\n",
    "s"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because we have a two terms relation,\n",
    "for a unique solution, we need to initial conditions,\n",
    "as an additional dictionary argument."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle - \\frac{\\sqrt{5} \\left(\\frac{1}{2} - \\frac{\\sqrt{5}}{2}\\right)^{n}}{5} + \\frac{\\sqrt{5} \\left(\\frac{1}{2} + \\frac{\\sqrt{5}}{2}\\right)^{n}}{5}$"
      ],
      "text/plain": [
       "-sqrt(5)*(1/2 - sqrt(5)/2)**n/5 + sqrt(5)*(1/2 + sqrt(5)/2)**n/5"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "s2 = rsolve(Fib, f(n), {f(0): 0, f(1): 1})\n",
    "s2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us now compute the first 10 Fibonacci numbers.  For this, we return to Sage and convert the SymPy object into a symbolic expression with ``SR()``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<class 'sympy.core.add.Add'>"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "type(s2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<class 'sage.symbolic.expression.Expression'>"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fs2 = SR(s2)\n",
    "type(fs2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can evaluate and expand the ``sqrt(5)`` expressions,\n",
    "which are then automatically simplified."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0, 1, 1, 2, 3, 5, 8, 13, 21, 34]"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[fs2(n=k).expand() for k in range(10)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4. Finding the Nearest Algebraic Number"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With ``PSLQ`` one can compute the nearest algebraic number to a given floating-point approximation.  This is available in ``mpmath``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.41421356237310"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "s2 = RR(sqrt(2))\n",
    "s2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "from mpmath import pslq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1.00000000000000, 1.41421356237310, 2.00000000000000]"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pwrs = [s2**k for k in range(3)]\n",
    "pwrs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "PSLQ using prec 113 and tol 1.81899e-12\n",
      "0/100:  Error:     0.06   Norm: 0\n",
      "FOUND relation at iter 1/100, error: 2.0e-16\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[2, 0, -1]"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cff = pslq(pwrs, verbose=True)\n",
    "cff"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In ``cff`` are the coefficients of the polynomial ``2 - x^2``, which is the smallest polynomial for which ``s2`` is an approximate root."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath 10.3",
   "language": "sage",
   "name": "sagemath"
  },
  "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.10.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}