Lecture 36: Symbolic Computation with sympy
===========================================

SymPy is pure Python package for symbolic computation.
The *pure* means that, unlike in SageMath, no modifications to the 
Python scripting language have been made.
SymPy is integrated in SageMath and at times we need to be aware of
how to use its functionality properly.
SageMath is licensed under the GNU GPL, 
whereas the license of SymPy is BSD.

sympy outside SageMath
----------------------

Outside SageMath, we can compute with SymPy as in the Sage cell server,
point your browser to the URL <http://live.sympy.org>.

A comparison of SymPy versus SageMath is given at
<https://github.com/sympy/sympy/wiki/SymPy-vs.-Sage>.

SymPy is one of the core packages of the SciPy stack
<http://www.scipy.org>.

sympy inside SageMath
---------------------

We explicitly applied sympy to define a generator for the terms
in a Taylor series.

If we use sympy in SageMath, then we may have to be explicit 
about the use of its functions.
For example, the ``var`` in SageMath is different 
from the ``var`` in sympy.

.. index:: var

::

   import sympy
   x = sympy.var('x')
   type(x)

We see that the type of ``x`` is ``sympy.core.symbol.Symbol.``

Observe the difference with the ``var`` that we are normally using in SageMath:

::

   y = var('y')
   type(y)

and the type is
``sage.symbolic.expression.Expression.``

If we use the ``sqrt()`` on a variable that is a sympy symbol, 
then we get a symbolic expression of SageMath on return.

::

   print(type(sqrt(x)))
   print(type(sqrt(y)))

In both cases we see the same type.
We can make a generator object on the ``sqrt(1-x).``  That works.

::

   tsqrtx = sympy.series(sqrt(1-x),x,0,n=None)
   print(tsqrtx)
   print([next(tsqrtx) for k in range(10)])

But it will not work on ``sqrt(1-y).``

::

   tsqrty = sympy.series(sqrt(1-y),y,0,n=None)
   print(tsqrty)
   next(tsqrty)
   next(tsqrty)

Although ``tsqrty`` is also some type of generator,
the second application of the generator triggers 
the ``StopIteration`` exception.

So for the generating function for the Catalan numbers, 
we must use a sympy symbol ``x`` and not a SageMath variable ``y.``

::

   g = (1 - sqrt(1-4*x))/(2*x)
   cg = sympy.series(g,x,0,n=None)
   print([next(cg) for k in range(11)])
   print(catalan_number(10))


Pattern Matching
----------------

Another feature better in SymPy than in SageMath
is pattern matching,
which can be useful when combined with substitution.
To run this code we can also open a Terminal session, select Misc, 
and select python.

::

   reset()
   from sympy import *
   x = Symbol('x')
   y = Wild('y')
   d = (10*x**3).match(y*x**3)
   d

On return is a dictionary that has ``y_`` as key
with ``10`` as corresponding value, which expresses
the match between ``y`` and ``10``
in the expressions ``10*x**3`` and ``y*x**3``.

Solving Recurrence Relations
----------------------------

SymPy can solve recurrence relations.
Consider the :index:`Fibonacci numbers`,
defined as :math:`f(n) = f(n-1) + f(n-2)`.
We import ``rsolve`` and declare ``f`` as a function of ``n``.

.. index:: rsolve, Function

::

   reset()
   from sympy import Function, rsolve
   from sympy.abc import n
   f = Function('f')

The equation :math:`f(n) = f(n-1) + f(n-2)`
is entered as ``f(n+2) - f(n+1) - f(n)``
and assigned to the variable ``Fib``.
This variable is the first argument of ``rsolve``
and we solve for ``f(n)``.

::

   Fib = f(n+2) - f(n+1) - f(n)
   print('solving', Fib)
   s = rsolve(Fib, f(n))
   s

The solution ``s`` has the following form:
``C0*(1/2 + sqrt(5)/2)**n + C1*(-sqrt(5)/2 + 1/2)**n``.
Because we have a two terms relation,
for a unique solution, we need to initial conditions,
as an additional dictionary argument.

::

   s2 = rsolve(Fib, f(n), {f(0): 0, f(1): 1})
   s2

and the solution ``s2`` equals
``sqrt(5)*(1/2 + sqrt(5)/2)**n/5 - sqrt(5)*(-sqrt(5)/2 + 1/2)**n/5``.
As a verification, let us apply the formula to compute the
first 10 Fibonacci numbers.  But first, we need to convert the SymPy ``s2``
into a SageMath expression:

.. index:: SR, symbolic ring

::

   fs2 = SR(s2)

Now we can evaluate and expand the expressions involving ``sqrt(5)``,
which are then automatically simplified.

::

   [fs2(n=k).expand() for k in range(10)]

The Nearest Algebraic Number
----------------------------

SymPy depends only on ``mpmath`` and ``mpmath`` exports the ``pslq``
function which allows to compute the nearest algebraic number,
given a floating-point approximation.  Let us illustrate this
with \ :math:`\sqrt{2}` which is a root of ``x^2 - 2``.

::

   from mpmath import pslq
   s2 = RR(sqrt(2))
   pwrs = [s2**k for k in range(3)]

In ``pwrs`` are the first three powers of ``s2`` starting at the zeroth power.

::

   cff = pslq(pwrs, verbose=True)
   cff 

The output in ``cff`` is the list ``[2, 0, -1]`` which are the
coefficients of the polynomial \ :math:`2 - x^2`.
This is the polynomial for which ``s2`` is an approximate roots
with 53 bits of precision.

Assignments
-----------

1. What happens if you ask SageMath or sympy to compute a Taylor series
   expansion of ``sqrt(x)`` at ``x = 0``?  Explain the result.

2. The cost :math:`T(n)` of a merge sort on a list of :math:`n` numbers is
   governed by the following recurrence relation:

   .. math::

      T(n) = 2 T\left(\frac{n}{2}\right) + n-1, \quad T(1) = 0.

   Use rsolve to find an explicit solution for :math:`T(n)`.
   Simplify so you can see :math:`T(n)` is :math:`O(n \log(n))`.