2.2.1   Symbols

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))

2.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. 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

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

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

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.)

2.2.3   Linear Algebra

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

2.2.4   Calculus

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.