GiNaC is Not a Computer Algebra System Chris Dams September 1, 2006 - - PowerPoint PPT Presentation

GiNaC is Not a Computer Algebra System

Chris Dams∗ September 1, 2006

∗Chris.Dams@mi.infn.it

A bit of history

Around 1998 the XLOOPs project had problems with their use the Maple CAS. Maple was not very suitable to support large projects. Major rewrites were necessary every new Maple version.

Use a standardized language

Use C++ and implement algebraic capabilities within it.

• No need to invent your own language.
• A mainstream language has good tools available.
• There is a standard library with good data structures.
• Support for modern programming techniques.
• Use it seamlessly with any other C/C++-library.
• Operator overloading.

GiNaCs classes

expairseq pseries relational matrix exprseq add ncmul function indexed basic ex lst tensor clifford mul idx color varidx basic ex wraps B A is from derived container class abstract class atomic class class template tensdelta power fderivative symbol constant numeric wildcard structure

... ...

Example code: integrating delta functions

The Dirac-delta function is a “function” that has the properties

• δ(x) is zero for x = 0;
• δ(x − x0) dx = 1;
• think of it as the mass density of a point particle.

For an integral containing a delta function we have

• f(x)δ(x − x0) dx = f(x0).

More generally

• f(x)δ(g(x)) dx =
• g(xi)=0

f(xi) |g′(xi)|.

Main program

int main() { symbol x("x"); ex e = integ(x,1)*delta(2*x-1)*exp(x); cout << "The expression " << e << endl; cout << "can be simplified into " << trydeltaintegrals(e) << endl; return 0; } The output is: The expression delta(-1+2*x)*exp(x)*integ(x,1) can be simplified into 1/2*exp(1/2)

The function that makes it possible

DECLARE_FUNCTION_2P(integ); REGISTER_FUNCTION(integ, dummy()); DECLARE_FUNCTION_1P(delta); REGISTER_FUNCTION(delta, dummy()); ex trydeltaintegrals(const ex &x) { if (is_a<mul>(x)) ... code for handling multiplications ... if (is_a<add>(x)) return x.map(trydeltaintegrals); return x; }

The code for handling multiplications

for (int i=0; i<x.nops(); ++i) { if (!is_ex_the_function(x.op(i), integ)) continue; ex integ_var = x.op(i).op(0); for (int j=0; j<x.nops(); ++j) { if (!is_ex_the_function(x.op(j), delta)) continue; ... more checks and returning result ... } } return x;

More checks and returning result

ex delta_arg = x.op(j).op(0); if (!delta_arg.is_polynomial(integ_var)) continue; if (delta_arg.degree(integ_var) != 1) continue; ex a = delta_arg.coeff(integ_var, 1); ex b = delta_arg.coeff(integ_var, 0); ex result = x/x.op(i)/x.op(j); result = result.subs( integ_var == -b/a ); result = result * power(abs(a), -x.op(i).op(1)); result = trydeltaintegrals(result); return result;

Features

• Arbitrary precision arithmetic for integers, fractions and reals;
• data structures: lst, exvector, exmap, exhashmap;
• linear algebra: Linear systems, matrix operations;
• algebraic substitution:

apbqarbs → ap+rbq+s by using the method .subs( power(wild(0),wild(1)) * power(wild(0),wild(2)) == power(wild(0), wild(1)+wild(2)), subs_options::algebraic );

• polynomial arithmetic: expanding, collecting, division, GCD,

LCM, resultant, square-free decomposition;

• analysis: Symbolic differentiation/Series expansion;
• HEF: Indexed objects, clifford algebra, lie algebras of SU(2) and

SU(3);

• automatic dummy index renaming:

f2.subs(f == x_i*y_i) → xiyixjyj;

• generic clifford algebras;
• taking real and imaginary parts of expressions;
• numeric integration;
• series expansion of symbolic integrals.

Fundamental concepts

• Reference counting for memory management;
• most copies are shallow;
• automatic evaluation:

– x + x → 2x; – ex::ex(const basic &other) for conversion; – calls virtual basic::eval() methods; – only operations that are of complexity N log(N)p;

ptr<basic> ex::construct_from_basic(const basic & other) { if (!(other.flags & status_flags::evaluated)) { const ex & tmpex = other.eval(1); if ((other.get_refcount() == 0) && (other.flags & status_flags::dynallocated)) delete &other; return tmpex.bp; } else { if (other.flags & status_flags::dynallocated) { return ptr<basic>(const_cast<basic &>(other)); } else { basic *bp = other.duplicate(); bp->setflag(status_flags::dynallocated); return bp; } } }

Fundamental Concept (continued)

• built-in methods use virtual functions (extendable!);
• map functions & type switch;
• hashing (Fibonacci hash) for fast comparison.

Projects using GiNaC

• Feelfem (finite element method (FEM) code generator);
• gTybalt (combines GiNaC with TeXmacs and Root into an

interactive computer algebra system);

• Nested sums (supports for some kind of transcendental

functions);

• PURRS (C++ library for solving recurrence relations);
• PyGiNaC (Python interface to GiNaC);
• Swiginac (Python interface to GiNaC);
• Antipodes (computes divergent terms arising in

HEF-computations);

• Ecco (simulation of (electro)chemical reaction kinetics);
• ORSA (Orbit Reconstruction, Simulation and Analysis);
• SyFi (finite element method).

Get it from

http://www.ginac.de Also included in most Linux distributions.

Authors/contributors

• Founding fathers: Christian Bauer, Alexander Frink,

Richard Kreckel;

• contributers: Roberto Bagnara, Jonathan Brandmeyer,

Matti Peltom¨ aki, Do Hoang Son, Bernard Parisse, Pearu Peterson, Ben Sapp, Stefan Weinzierl, Oliver Welzel

• active developers: Christian Bauer, Chris Dams, Vladimir Kisil,

Richard Kreckel, Alexei Sheplyakov, Jens Vollinga.