******************************************
Symbolic Computation with GiNaC and Python
******************************************

.. Bibliographic Fields

:Version:   1.3.5
:Status:    draft
:Date:      2007-02-01
:Copyright: 2006 Guenter Milde.
            Released under the terms of the GNU General Public License 
            (v. 2 or later)
:Abstract:  Using the GiNaC open framework for symbolic computation within 
            the Python programming language

.. :Author: Guenter Milde
.. :Organization: organization.
.. :Contact: contact.
.. :Address: address.
.. :Dedication: topic.

.. sectnum::

.. contents::

Introduction
============

This tutorial is intended for the novice user who is new to GiNaC_ but
already has some background in Python programming. It is an adaption of
the `GiNaC tutorial`_ for doing Computer algebra with GiNaC_ and the
Python_ programming language using the swiginac_ wrapper package. 

This tutorial only covers the basics, the original `GiNaC tutorial`_
says
   
   ... since a hand-made documentation like this one is difficult to
   keep in sync with the development, the actual documentation is
   inside the sources in the form of comments. That documentation may
   be parsed by one of the many Javadoc-like documentation systems. If
   you fail at generating it you may access it from the GiNaC_ home
   page. It is an invaluable resource not only for the advanced user
   who wishes to extend the system (or chase bugs) but for everybody
   who wants to comprehend the inner workings of GiNaC. This little
   tutorial on the other hand only covers the basic things that are
   unlikely to change in the near future.

The examples in this tutorial assume a working installation of GiNaC
and the swiginac package (see the Installation_ chapter for how to get
there).

.. _`GiNaC tutorial`: http://www.ginac.de/tutorial/
.. _GiNaC: http://www.ginac.de
.. _Python: http://www.python.org
.. _swiginac: http://swiginac.berlios.de

Source
------

The `source of this tutorial`_ is written in `reStructured Text`_. To
translate it, you need the Python docutils_. With docutils installed,
the command ``rst2html ginac-tutorial.py.txt`` will produce the html
file, while ``rst2pdf.py ginac-tutorial.py.txt`` will produce a PDF
version.

The bidirectional text <-> code converter PyLit_ can be used to
transform the tutorial source to a Python module providing the example
functions and to run the examples in a doctest. (This is why it has the
``.py.txt`` double extension.)

.. _source of this tutorial:
  http://svn.berlios.de/svnroot/repos/swiginac/trunk/doc/ginac-tutorial.py.txt
.. _docutils: http://docutils.sourceforge.net/rst.html
.. _reStructured Text: http://docutils.sourceforge.net/rst.html
.. _pylit: http://pylit.berlios.de

A Tour of GiNaC
===============

This quick tour of GiNaC wants to arise your interest in the subsequent
chapters by showing off a bit.

Please excuse us if it leaves many open questions.

How to use it from within Python
--------------------------------

Neither GiNaC nor swiGiNaC define a language of their own as conventional
CAS do. Instead, swiGiNaC converts GiNaC into an open framework for symbolic
computation within the Python programming language

Here is how to generate and print a simple (and rather pointless) bivariate
polynomial with some large coefficients::

  from swiginac import *

  def bipol():
      x = symbol('x')
      y = symbol('y')
      
      poly = 0
      for i in xrange(3):
           poly += factorial(i+16) * x**i * y**(2-i)
      
      return poly

Now we can evaluate the definition and try it:

  >>> execfile("ginac-tutorial.py")
  >>> print bipol()
  355687428096000*y*x+20922789888000*y**2+6402373705728000*x**2

Next, there is a more meaningful function that generates Hermite polynomials
in a specified free variable::

  def HermitePoly(x, n):
       HKer = exp(-x**2)
       # uses the identity H_n(x) == (-1)^n exp(x^2) (d/dx)^n exp(-x^2)
       return normal((-1)**n * 1/HKer * diff(HKer, x, n))
   
When run, this will type out the Hermite polynomials as symbolic terms:

  >>> z = symbol("z");
  >>> for i in xrange(6):
  ...     print "H_%d(z) == %s"%(i, HermitePoly(z,i))
  H_0(z) == 1
  H_1(z) == 2*z
  H_2(z) == -2+4*z**2
  H_3(z) == -12*z+8*z**3
  H_4(z) == 12+16*z**4-48*z**2
  H_5(z) == 120*z+32*z**5-160*z**3

This method of generating the coefficients is of course far from optimal for
production purposes.


In order to show some more examples of what GiNaC can do, the original
tutorial uses the ``ginsh``, a simple GiNaC interactive shell. We will show
that the swiginac_ Python interface provides a similar convenient window
into GiNaC's capabilities with all of the power of Python_ added.


What it can do for you
----------------------


After starting a python interpreter, import swiginac with:

  >>> from swiginac import *

Now you can test and experiment with GiNaC's features much like in other
Computer Algebra Systems. 

While the simple ``ginsh`` shell does not provide programming constructs
like loops or conditionals, ``swiginac`` combines the symbolic capacities of
GiNaC with the power of the full grown and established Python programming
language.

One main difference between "normal" CAS systems (or ginsh) and computer
algebra with Python is the need to define any symbols before use. Otherwise
Python will complain about unidentified variables:

  >>> sin(x)
  Traceback (most recent call last):
    File "<stdin>", line 1, in ?
  NameError: name 'x' is not defined

  >>> x = symbol('x')
  >>> sin(x)
  sin(x)

OTOH, it is easy to use the power of Python to define a set of
symbols, e.g. in a script::

  # define all lowercase letters (a-z) as Symbols
  import sys, string

  def new_symbol(*names):
      """Create new symbolic variables for all names in the argument list"""
      for name in names:
          setattr(sys.modules[__name__], name, symbol(name))

and use this later as e.g.:

  >>> for name in string.lowercase:
  ...    new_symbol(name)

  >>> type(u)
  <class 'swiginac.symbol'>

The ``numeric`` class can manipulate arbitrary precision integers in a very
fast way. Rational numbers are automatically converted to fractions of
coprime integers. As by default Python will convert a number to one of its
internal numerical data types, it is necessary to specify the wish to treat
a number as arbitrary precision integer. In an expression, it is often
sufficient to define just one numeric value as instance of the numeric class:

  >>> x = numeric(3)**150
  >>> print x
  369988485035126972924700782451696644186473100389722973815184405301748249

  >>> y = numeric(3)**149
  >>> print y
  123329495011708990974900260817232214728824366796574324605061468433916083

  >>> print x/y
  3

  >>> print y/x
  1/3

When specifying rational numbers, the expression that is the argument
to `numeric` is evaluated using Pythons "standard" arithmetic, which is
probabely not what we want:

  >>> print numeric(1/3)
  0

  >>> print numeric(1./3)
  0.33333333333333331483

To prevent evaluation by the "normal" Python operators, rational numbers can
be input as string value or as (numerator, denominator) tuple:

  >>> print numeric("1/3")
  1/3

  >>> print numeric(1, 3)
  1/3

  >>> print numeric(3, 9)
  1/3

Exact numbers are always retained as exact numbers and only evaluated as
floating point numbers if requested. For instance, products of sums of
radicals can be expanded:

  >>> a = symbol('a')
  >>> expand((1+a**numeric('1/5')-a**numeric('2/5'))**3)
  1-a**(6/5)+3*a**(1/5)-5*a**(3/5)+3*a

Numeric radicals are dealt pretty much as symbols:

  >>> expand((1+3**numeric('1/5')-3**numeric('2/5'))**3)
  10-5*3**(3/5)

  >>> ((1+3**numeric('1/5')-3**numeric('2/5'))**3).evalf()
  0.33408977534118624228


The ``evalf()`` method used above converts any number in GiNaC's expressions
into floating point numbers:

  >>> numeric('1/7').evalf()
  0.14285714285714285714

.. note:: In GiNaC or ginsh, this can be done to arbitrary predefined
    accuracy:

   > Digits=150;
   150
   > evalf(1/7);
   0.1428571428571428571428571428571428571428571428571428571428571428571428

 I did not find an equivalent setting in swiginac, though. Either this is
 not available (yet) or I just did not find it.

The return value of ``evalf()`` is still a ``numeric`` value:

  >>> type(numeric('1/7').evalf())
  <class 'swiginac.numeric'>

You can convert it to a Python datatype by the usual Python means or 
using the ``to_double()``, ``to_int()``, or ``to_long()`` methods:

  >>> type(numeric('1/7').to_double())
  <type 'float'>


Exact numbers that can be manipulated in GiNaC include predefined constants
like Archimedes' Pi. They can both be used in symbolic manipulations (as an
exact number) as well as in numeric expressions (as an inexact number):

  >>> x = symbol('x')
  >>> a = Pi**2 + x
  >>> print a
  Pi**2+x

  >>> a.evalf()
  9.869604401089358619+x

  >>> a.subs(x == 2).evalf()
  11.869604401089358619


.. note:: The ``subs()`` method is used, as changing the value of ``x`` will
          not propagate to the already defined expression ``a``:

               >>> x = 2
               >>> a.evalf()
               9.869604401089358619+x

Built-in functions evaluate immediately to exact numbers if this is
possible. Conversions that can be safely performed are done immediately;
conversions that are not generally valid are not done:

  >>> x = symbol('x')
  >>> cos(42*Pi)
  1
  >>> cos(acos(x))
  x
  >>> acos(cos(x));
  acos(cos(x))

(Note that converting the last input to x would allow one to conclude that
42*Pi is equal to 0.)


Linear equation systems can be solved along with basic linear algebra
manipulations over symbolic expressions. swiginac offers a matrix
class for this purpose:

  
  >>> x = symbol('x')
  >>> y = symbol('y')
  >>> z = symbol('z')
  >>> a = symbol('a')
  >>> b = symbol('b')

  >>> lsolve(a+x*y==z,x)
  y**(-1)*(z-a)

  >>> lsolve([3*x+5*y == 7, -2*x+10*y == -5], [x, y])
  [x==19/8, y==-1/40]

  >>> M = matrix([[1, 3], [-3, 2]])
  >>> determinant(M);
  11

  >>> l = symbol('lambda', '\lambda')  # lambda is a Python keyword
  >>> charpoly(M,l)
  11-3*lambda+lambda**2
  
  >>> A = matrix([[1, 1], [2, -1]])
  >>> A + 2*M
  [[1,1],[2,-1]]+2*[[1,3],[-3,2]]
  >>> evalm(A + 2*M)
  [[3,7],[-4,3]]

  >>> B = matrix([[0, 0, a], [b, 1, -b], [-1/a, 0, 0]])
  >>> evalm(B**(2**12))
  [[1,0,0],[0,1,0],[0,0,1]]

Multivariate polynomials and rational functions may be expanded, collected
and normalized (i.e. converted to a ratio of two coprime polynomials):

  >>> a = x**4 + 2*x**2*y**2 + 4*x**3*y + 12*x*y**3 - 3*y**4
  >>> b = x**2 + 4*x*y - y**2

  >>> expand(a*b)
  3*y**6+8*y*x**5+16*y**3*x**3+x**6+43*y**4*x**2-24*y**5*x+17*y**2*x**4
  >>> collect(a+b,x)
  (1+2*y**2)*x**2-3*y**4+x**4+x*(4*y+12*y**3)+4*y*x**3-y**2
  >>> collect(a+b,y)
  12*y**3*x-3*y**4+x**4+(-1+2*x**2)*y**2+x**2+y*(4*x+4*x**3)

  >>> normal(a/b)
  x**2+3*y**2

You can differentiate functions and expand them as Taylor or Laurent series
in a very natural syntax (the second argument of series is a relation
defining the evaluation point, the third specifies the order):

  >>> diff(tan(x),x)
  1+tan(x)**2

  >>> series(sin(x),x==0,4)
  1*x+(-1/6)*x**3+Order(x**4)

  >>> series(1/tan(x),x==0,2)
  1*x**(-1)+(-1/3)*x+Order(x**2)

  series(tgamma(x),x==0,2)
  1*x**(-1)+(-Euler)+(1/12*Pi**2+1/2*Euler**2)*x+Order(x**2)

  series(tgamma(x),x==0,2).evalf()
  1.0*x**(-1)+(-0.5772156649015328606)+0.9890559953279725555*x+Order(x**2)

  >>> series(tgamma(2*sin(x)-2),x==Pi/2,3)
  (-1)*(-1/2*Pi+x)**(-2)+(-1/12-Euler)+(-1/240-1/12*Pi**2-1/2*Euler**2)*(-1/2*Pi+x)**2+Order((-1/2*Pi+x)**3)

Often, functions don't have roots in closed form. Nevertheless, GiNaC can
compute a solution numerically:

   TODO: is there an equivalent to fsolve in swiginac?


Symbolic types can always be used as tags for different types of objects.
Converting from "wrong" units to the metric system is an example:

  >>> new_symbol('m', 'kg', 'lb', 'inch')
  >>> inch = .0254*m
  >>> lb = .45359237*kg
  >>> 200*lb/inch**2
  (140613.91592783187407)*kg*m**(-2)

BTW, the *Scientific.Physics.PhysicalQuantities* module from the Scientific_
package provides an elaborated environment specifically designed for working
with units.

.. _Scientific: http://dirac.cnrs-orleans.fr/ScientificPython/


Installation
============

Installation is covered in the original documentation of the relevant
packages. Instead of repeating it, just pointers will be given.

GiNaC
-----

The GiNaC framework for symbolic computation within the C++ programming
language (and its tutorial) is Copyright 1999-2006 Johannes Gutenberg
University Mainz, Germany. It is released under the 
`GNU General Public License`__ (GPL).

__ http://www.gnu.org/copyleft/gpl.html

GiNaC_ sources can be fetched from its home site under the section
`Download GiNaC now`_ which also says:

  Since GiNaC is packaged by a couple of software distributors, you may want
  to try a conveniently packaged binary first (but please don't send us
  bug-reports about these packages since they may be outdated). We are
  currently aware of packages in Debian__, SuSE, Mandrake, Fedora, and
  FreeBSD.
  
The installation__ is covered in the `GiNaC tutorial`_ (and in the INSTALL
file of the tarball).

.. _download GiNaC now:  http://www.ginac.de/Download.html
__ http://www.ginac.de/tutorial/Installation.html#Installation
__ http://packages.qa.debian.org/g/ginac.html


swiginac
--------

swiginac_ is a Python interface to GiNaC, built with SWIG_. The aim of
swiginac is to make the functionality of GiNaC accessible from Python as
an extension module.

Current status is beta; a lot (but not all yet) of the GiNaC classes are
exposed, virtually all of the GiNaC tests pass.

For more information, documentation and software downloads, visit the
`swiginac group pages on BerliOS`__.

The `source tarball`__ contains installation instructions in a
INSTALL.txt file.

.. _SWIG: http://www.swig.org/
__ http://developer.berlios.de/projects/swiginac/
__ http://prdownload.berlios.de/swiginac/swiginac-0.9.4.tgz


PyGiNaC
-------

PyGiNaC_ is another Python package that provides an interface to the C++
library GiNaC. PyGiNaC is implemented with the help of the `Boost.Python`_
library. At the moment, the package is more or less in an "alpha" state,
i.e. the GiNaC classes are only partially exposed, parts of regression test
suite are unconverted and no actual documentation exists. However, all of
this is supposed to change in the future. 


.. _PyGiNaC: http://pyginac.sourceforge.net/

.. _`Boost.Python`: http://www.boost.org/libs/python/doc/index.html

At the moment there are no released tarballs of PyGiNaC. Instead, it can be
downloaded directly from the CVS repository at SourceForge. (See the
PyGiNaC_ page for instructions.)
