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 with the Python programming language
Contents
This tutorial is intended for the user who is new to GiNaC but already has some background in Python programming. It is an adaption of the GiNaC tutorial (for symbolic computation within the C++ programming language) for 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).
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.
GiNaC does not define a programming language of is own as conventional computer algebra systems (CAS) do. GiNaC is a C++ library for applications in need of symbolic manipulation.
Python is such an application.
Swiginac is a Python interface for GiNaC adding symbolic mathematics to the many features of Python. After starting a python interpreter, import the swiginac module with:
>>> from swiginac import *
Now you can test and experiment with GiNaC's features much like in other Computer Algebra Systems.
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 import the function and try it out:
>>> from swiginac_tutorial import bipol >>> 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:
>>> from swiginac_tutorial import HermitePoly >>> 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 GiNaC at work, the original tutorial uses the GiNaC interactive shell ginsh. The simple ginsh shell does not provide programming constructs like loops or conditionals, whereas swiginac combines the symbolic capacities of GiNaC with the ease of use and the power of the full grown and established programming language Python.
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)
You can also define a set of symbols in a loop:
>>> import sys >>> for sym in ['x', 'y', 'z']: ... setattr(sys.modules[__name__], sym, symbol(sym))
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+3*a**(1/5)-5*a**(3/5)
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
There is no equivalent setting in swiginac (yet).
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
>>> a x+Pi**2 >>> a.evalf() 9.869604401089358619+x >>> a.subs(x == 2).evalf() 11.869604401089358619
Note
In the last example, 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'); l = symbol('l')
>>> lsolve(a+x*y==z,x) (z-a)*y**(-1) >>> 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 >>> charpoly(M,l) 11+l**2-3*l
>>> A = matrix([[1, 1], [2, -1]]) >>> 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) 43*y**4*x**2+3*y**6-24*y**5*x+17*y**2*x**4+8*y*x**5+x**6+16*y**3*x**3 >>> collect(a+b,x) -3*y**4+x*(4*y+12*y**3)+(1+2*y**2)*x**2+4*y*x**3+x**4-y**2 >>> collect(a+b,y) -3*y**4+x**2+(4*x**3+4*x)*y+12*y**3*x+x**4+y**2*(-1+2*x**2)
>>> 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)*(x-1/2*Pi)**(-2)+(-1/12-Euler)+(-1/240-1/12*Pi**2-1/2*Euler**2)*(x-1/2*Pi)**2+Order((x-1/2*Pi)**3)
Often, functions don't have roots in closed form. Nevertheless, it's quite easy to compute a solution numerically:
>>> fsolve(cos(x)==x,x,0,2) 0.73908513321516064166 >>> f=exp(sin(x))-x >>> X=fsolve(f,x,-10,10) >>> X 2.2191071489137460327 >>> subs(f,x==X) -2.168404344971008868E-19
Notice how the final result above differs slightly from zero by about 2*10^(-19). This is because with 20 decimal digits precision the root cannot be represented more accurately than X. Such inaccuracies are to be expected when computing with finite floating point values.
Symbolic types can always be used as tags for different types of objects. Converting from "wrong" units to the metric system is an example:
>>> m = symbol('m'); kg = symbol('kg') >>> inch = .0254*m >>> lb = .45359237*kg >>> 200*lb/inch**2 (140613.91592783187407)*m**(-2)*kg
However, as generic symbols do not know about special properties of physical units, it can not simplify all expressions with unit symbols:
>>> print abs(-3), abs(-3*m) 3 abs(-3*m)
(We know, that m is positiv definit, but GiNaC does not.)
The Scientific.Physics.PhysicalQuantities module from the Scientific package provides an elaborated environment specifically designed for working with units. Unfortunately, it does not collaborate well with swiginac.
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).
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 tarballs contain installation instructions in an INSTALL.txt file.
The source of this tutorial is written in reStructured Text. To translate it, you need the Python docutils. With docutils installed, the command rst2html swiginac_tutorial.py.txt will produce the html file, while rst2pdf.py swiginac_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.)