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

2.1 How to use it from within C++

The GiNaC open framework for symbolic computation within the C++ programming language does not try to define a language of its own as conventional CAS do. Instead, it extends the capabilities of C++ by symbolic manipulations. Here is how to generate and print a simple (and rather pointless) bivariate polynomial with some large coefficients:

#include <iostream>
#include <ginac/ginac.h>
using namespace std;
using namespace GiNaC;

int main()
{
    symbol x("x"), y("y");
    ex poly;

    for (int i=0; i<3; ++i)
        poly += factorial(i+16)*pow(x,i)*pow(y,2-i);

    cout << poly << endl;
    return 0;
}

Assuming the file is called `hello.cc', on our system we can compile and run it like this:

$ c++ hello.cc -o hello -lcln -lginac
$ ./hello
355687428096000*x*y+20922789888000*y^2+6402373705728000*x^2

(See section B. Package Tools, for tools that help you when creating a software package that uses GiNaC.)

Next, there is a more meaningful C++ program that calls a function which generates Hermite polynomials in a specified free variable.

#include <iostream>
#include <ginac/ginac.h>
using namespace std;
using namespace GiNaC;

ex HermitePoly(const symbol & x, int n)
{
    ex HKer=exp(-pow(x, 2));
    // uses the identity H_n(x) == (-1)^n exp(x^2) (d/dx)^n exp(-x^2)
    return normal(pow(-1, n) * diff(HKer, x, n) / HKer);
}

int main()
{
    symbol z("z");

    for (int i=0; i<6; ++i)
        cout << "H_" << i << "(z) == " << HermitePoly(z,i) << endl;

    return 0;
}

When run, this will type out

H_0(z) == 1
H_1(z) == 2*z
H_2(z) == 4*z^2-2
H_3(z) == -12*z+8*z^3
H_4(z) == -48*z^2+16*z^4+12
H_5(z) == 120*z-160*z^3+32*z^5

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 we will now use the ginsh, a simple GiNaC interactive shell that provides a convenient window into GiNaC's capabilities.

2.2 What it can do for you

After invoking ginsh one can test and experiment with GiNaC's features much like in other Computer Algebra Systems except that it does not provide programming constructs like loops or conditionals. For a concise description of the ginsh syntax we refer to its accompanied man page. Suffice to say that assignments and comparisons in ginsh are written as they are in C, i.e. = assigns and == compares.

It can manipulate arbitrary precision integers in a very fast way. Rational numbers are automatically converted to fractions of coprime integers:

> x=3^150;
369988485035126972924700782451696644186473100389722973815184405301748249
> y=3^149;
123329495011708990974900260817232214728824366796574324605061468433916083
> x/y;
3
> y/x;
1/3

Exact numbers are always retained as exact numbers and only evaluated as floating point numbers if requested. For instance, with numeric radicals is dealt pretty much as with symbols. Products of sums of them can be expanded:

> expand((1+a^(1/5)-a^(2/5))^3);
1+3*a+3*a^(1/5)-5*a^(3/5)-a^(6/5)
> expand((1+3^(1/5)-3^(2/5))^3);
10-5*3^(3/5)
> evalf((1+3^(1/5)-3^(2/5))^3);
0.33408977534118624228

The function evalf that was used above converts any number in GiNaC's expressions into floating point numbers. This can be done to arbitrary predefined accuracy:

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

Exact numbers other than rationals 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):

> a=Pi^2+x;
x+Pi^2
> evalf(a);
9.869604401089358619+x
> x=2;
2
> evalf(a);
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:

> 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. In C++ GiNaC offers a matrix class for this purpose but we can see what it can do using ginsh's bracket notation to type them in:

> 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 = [ [1, 3], [-3, 2] ];
[[1,3],[-3,2]]
> determinant(M);
11
> charpoly(M,lambda);
lambda^2-3*lambda+11
> A = [ [1, 1], [2, -1] ];
[[1,1],[2,-1]]
> A+2*M;
[[1,1],[2,-1]]+2*[[1,3],[-3,2]]
> evalm(%);
[[3,7],[-4,3]]
> B = [ [0, 0, a], [b, 1, -b], [-1/a, 0, 0] ];
> evalm(B^(2^12345));
[[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;
12*x*y^3+2*x^2*y^2+4*x^3*y-3*y^4+x^4
> b = x^2 + 4*x*y - y^2;
4*x*y-y^2+x^2
> expand(a*b);
8*x^5*y+17*x^4*y^2+43*x^2*y^4-24*x*y^5+16*x^3*y^3+3*y^6+x^6
> collect(a+b,x);
4*x^3*y-y^2-3*y^4+(12*y^3+4*y)*x+x^4+x^2*(1+2*y^2)
> collect(a+b,y);
12*x*y^3-3*y^4+(-1+2*x^2)*y^2+(4*x+4*x^3)*y+x^2+x^4
> normal(a/b);
3*y^2+x^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);
tan(x)^2+1
> series(sin(x),x==0,4);
x-1/6*x^3+Order(x^4)
> series(1/tan(x),x==0,4);
x^(-1)-1/3*x+Order(x^2)
> series(tgamma(x),x==0,3);
x^(-1)-Euler+(1/12*Pi^2+1/2*Euler^2)*x+
(-1/3*zeta(3)-1/12*Pi^2*Euler-1/6*Euler^3)*x^2+Order(x^3)
> evalf(%);
x^(-1)-0.5772156649015328606+(0.9890559953279725555)*x
-(0.90747907608088628905)*x^2+Order(x^3)
> series(tgamma(2*sin(x)-2),x==Pi/2,6);
-(x-1/2*Pi)^(-2)+(-1/12*Pi^2-1/2*Euler^2-1/240)*(x-1/2*Pi)^2
-Euler-1/12+Order((x-1/2*Pi)^3)

Here we have made use of the ginsh-command % to pop the previously evaluated element from ginsh's internal stack.

If you ever wanted to convert units in C or C++ and found this is cumbersome, here is the solution. Symbolic types can always be used as tags for different types of objects. Converting from wrong units to the metric system is now easy:

> in=.0254*m;
0.0254*m
> lb=.45359237*kg;
0.45359237*kg
> 200*lb/in^2;
140613.91592783185568*kg*m^(-2)