SwiGiNaC: Symbolic Computation with Python

Basic Examples

Contents

1   Import

This script will give a basic overview of what is possible combining symbolic computation by GiNac with the flexibility and richness of Python

To load swiginac, use one of the import variants:

import swiginac
from swiginac import *

Objects defined in this tutorial can be imported with

>>> from swiginac_basics import *

The swiginac module wraps functions, data and classes defined in the ginac C++ library in a python module. There are more than 600 objects

>>> len(dir(swiginac))
661

but only the most basic will be described in this tutorial.

Completing the swiginac functions and classes with docstrings is an open task. Some introspection is possible with e.g.

# from pprint import pprint
# pprint(dir(swiginac))
# pprint(vars(swiginac))

2   Objects

2.1   Symbols

A symbolic indeterminante or symbolic variable is a placeholder for a value in an expression.

Symbols are basic units in Swiginac:

>>> a = symbol('a')
>>> a
a

The datatype is:

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

It defines two methods but inherits more than 100 others:

>>> len(dir(a))
108

If you want to list them all, try pprint(dir(a))

In most Computer Algebra Systems (CAS), any unbound variable can be used as a symbol. In Python, you must define a variable as symbol before you can use it in expressions.

>>> del(x)
>>> y = 3*x
Traceback (most recent call last):
  File "<stdin>", line 1, in ?
NameError: name 'x' is not defined
>>> x = symbol('x')
>>> y = 3*x

y is now an expression:

>>> y
3*x

In order to re-use y as a symbol, an ordinary CAS would require you to delete it. In Python, you overwrite its current binding with a new assignment to a symbol instance:

>>> y = symbol('y')
>>> y
y

2.1.1   Defining a set of symbols

If we initialize a set of symbols, we might want to use a loop rather than a lot of lines with individual assignments.

We can use the feature that the dictionary returned by the built-in function globals can be manipulated

>>> for name in ['gamma', 'delta', 'epsilon']:
...     globals()[name] = swiginac.symbol(name)

To define the small latin letters a-z as symbols in this module,

import string as _string

for name in _string.lowercase:
  globals()[name] = swiginac.symbol(name)

Which results in

>>> print type(delta), type(x)
<class 'swiginac.symbol'> <class 'swiginac.symbol'>

2.2   Numbers

The numeric class can manipulate arbitrary precision integers in a very fast way.

Rational numbers are automatically converted to fractions of coprime integers.

>>> numeric('4/12')
1/3

The expression given to numeric is evaluated using Pythons "standard" arithmetic, which is probabely not what we want when specifying rational numbers:

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

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

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

Often, it is sufficient to specify one number in an expression as numeric:

>>> numeric(1)/3
1/3

2.2.1   Converting numbers

The evalf method converts any number in GiNaC's expressions into floating point numbers:

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

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 using the to_double, to_int, or to_long methods:

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

There are other converting methods, like numeric.evalm or numeric.to_cl_N but what do they do?

2.2.2   Complex Numbers

Complex numbers are expressad with the imaginary unit I

>>> I, I**2, 2+3*I
(I, -1, 2+3*I)

However, they are numbers (even if they look like a symbol, sum, or product):

>>> type(I), type(2 + 3*I)
(<class 'swiginac.numeric'>, <class 'swiginac.numeric'>)

Python's complex numbers can currently not be converted to numeric objects:

>>> z_gi = numeric(2+3j)
Traceback (most recent call last):
    ...
NotImplementedError: Wrong number of arguments for overloaded function 'new_numeric'.
  Possible C/C++ prototypes are:
    GiNaC::numeric(int)
    GiNaC::numeric(unsigned int)
    GiNaC::numeric(long)
    GiNaC::numeric(unsigned long)
    GiNaC::numeric(long,long)
    GiNaC::numeric(double)
    GiNaC::numeric(char const *)
    GiNaC::numeric(cln::cl_N const &)
<BLANKLINE>

A workaround is to use an expression:

z_py = 2+3j
z_gi = z_py.real + z_py.imag*I

which will be simplified to a number:

>>> print z_gi, type(z_gi)
2.0+3.0*I <class 'swiginac.numeric'>

How do complex expression evaluate?

evalf converts real and imaginary part to floating point numbers

>>> z_gi.evalf()
2.0+3.0*I

to_double() returns the real part as double instance

>>> z_gi.to_double()
2.0

While Phyton's complex class stores real and imaginary part as attributes, numeric provides methods to retrieve them (as numeric instances):

>>> z_py.real, z_gi.real
(2.0, <bound method numeric.real of 2.0+3.0*I>)
>>> z_gi.real(), z_gi.imag()
(2.0, 3.0)

2.3   Constants

Some mathematical constants are available as well

>>> print Pi, type(Pi), Pi.evalf()
Pi <class 'swiginac.constant'> 3.1415926535897932385
>>> print Euler, type(Euler), Euler.evalf()
Euler <class 'swiginac.constant'> 0.5772156649015328606

Swiginac functions know some simplifications for special values

>>> print sin(Pi), cos(Pi)
0 -1

2.4   Expressions

Expressions in GiNaC are built with symbols and numbers.

They are converted to a "canonical" representation and have a class that depends on the expression:

ex1 = a/2 + b
ex2 = (a+b)/2
ex3 = (a+b)/c
>>> print ex1, type(ex1)
1/2*a+b <class 'swiginac.add'>
>>> print ex2, type(ex2)
1/2*a+1/2*b <class 'swiginac.add'>
>>> print ex3, type(ex3)
(a+b)*c**(-1) <class 'swiginac.mul'>

In the Symbolic package, there is the common class Symbolic.Expr with some additional methods.

2.5   Functions

About 40 mathematical functions are available in swiginac. All trigonometric and hyperbolic functions are implemented. For a full list of available functions see the file doc/examples/functions.py.

Some functions are simplified if this is mathemetically sound and non-ambiguous.

>>> sin(x), sin(0), sin(Pi/2)
(sin(x), 0, 1)
>>> cos(x), cos(0), cos(Pi/2)
(cos(x), 1, 0)
>>> tan(x), tan(0), tan(Pi/2)
Traceback (most recent call last):
  ...
StandardError: tan_eval(): simple pole

But not all simplifications are implemented

>>> sin(x)**2 + cos(x)**2
sin(x)**2+cos(x)**2
z_s = cos(x) + I*sin(x)
z_e = exp(I*x)
>>> z_s/z_e
(cos(x)+I*sin(x))*exp(I*x)**(-1)

Is there a way get more simplifications?

2.6   Matrices

Marices can be defined giving the number of rows and columns:

M1 = matrix(2,2)

Elements will be initialized with 0 and can be set individually:

M1[0,0] = x
M1[1,1] = y
>>> M1
[[x,0],[0,y]]
>>> M1.printc()
'[[x,0.0],[0.0,y]]'
>>> M1.printlatex()
'\\left(\\begin{array}{cc}x&0\\\\0&y\\end{array}\\right)'

Alternatively, they can be initialized by a list of lists:

M2 = matrix([[x,0],[0, y]])

Diagonal matrices can be created with a special constructor:

M3 = diag_matrix([x, y])
>>> bool(M1 == M2 == M3)
True

Vectors are (1,n) or (n,1) matrices:

A = matrix([range(4)])
B = matrix([[a], [b], [c]])
C = matrix([[0], [1], [2]])

But why is the matrix product not working?

>>> print A*B
[[0,1,2,3]]*[[a],[b],[c]]
>>> C = matrix([[0], [1], [2]])
>>> print A*C
[[0,1,2,3]]*[[0],[1],[2]]

Functions are not broadcast to the matrix elements as in NumPy

>>> sin(A)
sin([[0,1,2,3]])

How about substitution?

>>> sin(x).subs(x == A)
sin([[0,1,2,3]])

3   Symbolic calculations

3.1   Differentiation

Objects have the method diff for differentiation which is also called by the diff function:

>>> P = x**5 + x**2 + x
>>> P.diff(x) == diff(P, x)
1+5*x**4+2*x==1+5*x**4+2*x
>>> P.diff(x, 2) == diff(P, x, 2)
2+20*x**3==2+20*x**3
>>> diff(sin(exp(x)), x) == sin(exp(x)).diff(x)
cos(exp(x))*exp(x)==cos(exp(x))*exp(x)

3.2   Simple integral support

We can construct integral objects and integrate either symbolically or numerically:

x = symbol('x')
integ = integral(x, 0, 1, x*x)
>>> integ
[integral object]
>>> print integ.printlatex()
\int_{0}^{1} dx\,x^{2}
>>> integ.eval_integ()
1/3
>>> integ.evalf()
0.33333333333333333332

3.3   Substitution

Algebraic objects in expressions can be substituted:

s0 = sin(exp(x))
s1 = s0.subs(exp(x)==sqrt(y))
s2 = s1.subs(y==2)
>>> s1
sin(y**(1/2))
>>> s2.evalf()
0.98776594599273552706
>>> float(s2.evalf())
0.98776594599273548

It is possible to replace individual symbols or terms. Several substitutions can be given in a sequence (list or tuple):

s3 = sin(x+y+z)

s4 = s3.subs([x==1, y==2])
s5 = s3.subs(x+y+z==6)
>>> print s4
sin(3+z)
>>> print s5
sin(6)

But sub-expressions do not match:

s6 = s3.subs(x+y==4)
>>> s6
sin(x+y+z)
>>> s3.subs(x+y==4)
sin(x+y+z)

Strange: Some substitutions work in the module but fail the doctest:

>>> print s3.subs([x==1, y==2])
sin(3+z)
Failed example:
print s3.subs([x==1, y==2])
Expected:
sin(3+z)
Got:
sin(x+y+z)
>>> print s3.subs([x==1, y==2, z==3])
sin(6)
Failed example:
print s3.subs([x==1, y==2, z==3])
Expected:
sin(6)
Got:
sin(3+x+y)
>>> s7 = s3.subs([x==1, y==2, z==3])
>>> s7
sin(6)
Failed example:
s7
Expected:
sin(6)
Got:
sin(3+x+y)

3.4   Solving linear systems

The function lsolve solves linear systems:

lgs = [3*x + 5*y == 2, 5*x+y == -3]

ev = lsolve(lgs, [x,y])
>>> ev
[x==-17/22, y==19/22]
>>> lgs[0].subs(ev)
2==2
>>> lgs[1].subs(ev)
-3==-3

3.5   Taylor series expansion

Expressions can expand themselves as a Taylor series:

x =symbol("x")
taylor = sin(x).series(x==0, 8)
>>> taylor
1*x+(-1/6)*x**3+1/120*x**5+(-1/5040)*x**7+Order(x**8)

4   Output Formatting

4.2   Format the output of numbers

The output format for numbers is sometimes too verbose for practical use, especially with floats as precision is govened by DIGITS (20):

>>> from swiginac import *
>>> x = Pi.evalf()
>>> print x
3.1415926535897932385

Python's built in function round() can be used to round any object that can be converted to a float to a given accuracy. The result is, however, not always as expected:

>>> round(x, 4)
3.1415999999999999

This is caused by the internal use of a number base different from 10.

4.3   String formatting

In Python, you can use string formatting to get a specific notation or accuracy:

The '%f' specifier can be used for swiginac.numeric instances that are convertible to floats:

>>> print "%f, %f, %f" % (x, sin(x), exp(x))
3.141593, -0.000000, 23.140693

It fails on objects that cannot be converted to floats:

>>> print "%f" % Pi
Traceback (most recent call last):
  ...
TypeError: float argument required

Save programming would require to e.g. test the datatype with

>>> x.is_real()
True

or use a try/catch clause.

The formatting options are described in the Python documentation on string formatting operations.

Here some examples:

Precision:

>>> print "%.10f" % x
3.1415926536

... is limited to DIGITS internally

>>> print "%.60f" % x
3.141592653589793115997963468544185161590576171875000000000000

Leading sign and zeros:

>>> print "%+.1f" % x
+3.1
>>> print "%05.1f" % x
003.1

Scientific notation

>>> print "%.3E  %.2e" % (x, x)
3.142E+00  3.14e+00
>>> print "%g  %g" % (x, x**40)
3.14159  7.69121e+19

4.4   Interaction with Python

4.4.1   Series and sums

A series could be defined as a Python sequence, e.g.

>>> x = symbol('x')
>>> [x/n for n in range(1, 10)]
[x, 1/2*x, 1/3*x, 1/4*x, 1/5*x, 1/6*x, 1/7*x, 1/8*x, 1/9*x]
>>> sinc5 = [sin(phi)/phi for phi in range(1,5)]; sinc5
[sin(1), 1/2*sin(2), 1/3*sin(3), 1/4*sin(4)]

Attention, the loop indices are now integers

>>> n, phi
(9, 4)

so it might be a good idea to use separate conventions for naming symbols and indices.

As the standard function sum is overloaded for swiginac classes, it can be used on sequences of symbols or expressions:

>>> sum((x, y, z))
y+x+z
>>> sum([x/n for n in range(1, 10)])
7129/2520*x

Compute the sum over the list defined earlier in this section:

>>> sum(sinc5)
1/2*sin(2)+1/3*sin(3)+sin(1)+1/4*sin(4)