GluCat: Generic library of universal Clifford algebra templates

SourceForge.net project page
GluCat project page hosted at:

GluCat is a library of template classes which model the universal Clifford algebras over the field of real numbers, with arbitrary dimension and arbitrary signature. GluCat implements a model of each Clifford algebra corresponding to each non-degenerate quadratic form up to a maximum number of dimensions.

GluCat classes were originally designed to be used as template parameters for other template libraries, such as Blitz++, deal.II, Matrix Template Library and POOMA. These template libraries expect a numeric class which implements "appropriate numeric semantics". To provide these semantics, the GluCat interface matches float or complex<> as much as possible, but the GluCat classes must be different from float or complex<>, because they model a different algebra.

Interface

The namespace is glucat::, which stands for Generic library of universal Clifford algebra templates.

The main template classes are index_set<LO,HI>, framed_multi<Scalar_T,LO,HI> and matrix_multi<Scalar_T,LO,HI>. Scalar_T is a class which represents a number field, such as float or double. LO, an Index_T <= 0 and HI, an Index_T >= 0, are the lower and upper limits for indices in index sets. Testing so far has been done with Index_T as int and Scalar_T as float, double and long double.

The notation for index_set<LO,HI> looks like: {-2,3,4}, a set of integers, with templated limits. 0 is never a member of the set. Each index set represents the non-commutative product of the generators corresponding to each index in strictly increasing order. This is usually called a strictly increasing multi-index in the literature.

eg. {-2,3,4} == {-2}*{3}*{4}

Input treats index sets as sets, so {-2,3,4} == {4,3,-2}. Output is in increasing order.

A framedmultivector, a multivector expressed with respect to a basis, "knows" which subalgebra contains it. It is the subalgebra defined by the frame. A frame here is an ordered set of basis vectors which generates a subalgebra.

The notation for framed_multi<> and matrix_multi<> is the framed notation: 1+2{-3}+4{-2}+5{-3,-2,2}, ie a sum of terms, where each term consists of a coordinate and an index_set.

framed_multi<float,LO,HI> and matrix_multi<float,LO,HI> represent elements of subalgebras of R(HI,-LO).

{-2}*{-2} == -1. {2}*{2} == 1.

The framed notation automatically takes care of the signature of the generators. You don't need to specify a signature to use framed_multi<> or matrix_multi<>, ie. the value of a multivector in this notation is unambiguous.

The framed notation, with negative indices, is not new. Mark Ashdown uses negative indices in GA Package for Maple V. G. P. Wene uses negative indices in his paper on infinite dimensional Clifford algebras. Leo Dorst defines various products in terms of index sets in his paper on "Honing geometric algebras". Finally, and most importantly, Arvind Raja's design and prototype code (the basis for GluCat) uses index sets internally.

framed_multi<> and matrix_multi<> provide the usual arithmetic operators +, -, * and /, as well as % for left contraction, & for Hestenes inner product and ^ for outer product. They also provide operator[] for subscripting using an index set, and operator() for grading.

framed_multi<> and matrix_multi<> also provide the usual functions, conj, real, imag, abs, norm, pow, exp, sqrt, log, cos, acos, cosh, acosh, sin, asin, sinh, asinh, tan, atan, tanh, atanh, plus multivector functions even, frame, inv, involute, reverse, pure, quad, scalar etc.

The aim is to get the notation as close as possible to that used by float or complex<>, while not straying too far from current notations for Clifford algebras.

Operation

One of the beauties of the Clifford algebras is that F(P,Q) contains subalgebras isomorphic to F(p,q) for p<=P and q<=Q.

It is therefore OK to add and multiply multivectors from different subalgebras. The result is, in general, contained in the subalgebra which spans the operands.

eg. (2{-1}+3{-2}) + (4{1,2}+5{1,4}) == 2{-1}+3{-2}+4{1,2}+5{1,4}

The idea in using the GluCat classes is to make the range LO...HI as large as you will need for your program. You can use the same class, say framed_multi<double,-4,4> to perform calculations in R(0,2) (quaternions), R(3,0) (3D geometric algebra), R(0,3) (3D Clifford algebra used in Clifford analysis) and R(1,3) (spacetime algebra).

Even better, you can just let LO and HI take default values, which are large enough for most purposes. The defaults are based on sizeof(set_value_t), where set_value_t is defined as unsigned long. On the Intel i386 architecture, unsigned long is 32 bits, DEFAULT_LO is -16 and DEFAULT_HI is 16. So to accommodate the subalgebras listed above, you can just use framed_multi<double>. This incurs minor overhead in size and speed when compared to framed_multi<double,-4,4>. On AMD 64 or Compaq Alpha, unsigned long is 64 bits, DEFAULT_LO is -32 and DEFAULT_HI is 32.

Performance

Conversion between framed_multi<> and matrix_multi<> is necessary in many circumstances. matrix_multi<> cannot operate directly on operands in different subalgebras, so such operations work by converting the operands to a common subalgebra by conversion to and from framed_multi<>.

Past a threshold number of dimensions, this conversion is implemented using a generalized Fast Fourier Transform.

Implementation

GluCat adapts Arvind Raja's prototype code, with permission, for the basic operations and adds /, inv, sqrt and the transcendental functions listed above.

The class index_set<LO,HI> is based on std::bitset<HI-LO+1>.

The class framed_multi<Scalar_T,LO,HI> is based on std::map<index_set<LO,HI>,Scalar_T> or (optionally) on __gnu_cxx::hash_map< const index_set<LO,HI>,Scalar_T, hash <LO,HI> > or std::tr1::unordered_map< const index_set<LO,HI>,Scalar_T, hash <LO,HI> > >. An object of class framed_multi<> represents a multivector as a sum of terms, where each term consists of a basis element and a coordinate. You can also think of this as a map from the set of index sets to the set of coordinate values. Where necessary, framed_multi<> implementation uses std::pair<index_set<LO,HI>,Scalar_T> to represent terms.

The class matrix_multi<> consists of a matrix, plus an index_set<> which represents the frame which defines the subalgebra. matrix_multi<> is currently based on the ublas:: templates provided by the Boost uBLAS Library. GluCat currently uses either sparse matrices in row compressed format or dense matrices.

The class matrix_multi<> uses the matrix generators of the minimum real representation of R(p,q) as per Porteous and Lounesto. A single instance of the class generator_table<> is used to construct and store the generators. This class, based on std::map<>, maps from signatures (p,q) to vectors of matrix generators.

Division and inverse currently use LU decomposition, plus iterative refinement, as per Higham. The class matrix_multi<> uses framed_multi<> for most operations other than +, -, *, /.

Conversion from framed_multi<> to matrix_multi<> uses the smallest subalgebra which contains the framed_multi<>.

The fold() and unfold() operations map an index set within a subalgebra to and from the isomorphic subalgebra defined by the corresponding contiguous index set. An index set defines a subalgebra by defining a frame, which is an ordered set of basis vectors. The indices in the index set, in order, define this ordered set. In turn, the ordered set of basis elements generates the subalgebra. Thus you can think of fold() and unfold() as renumbering the basis vectors within a frame.

Eg. {-4,-2,1,3}.fold({-4,-2,1,3}) == {-2,-1,1,2}
{-2,3}.fold({-4,-2,1,3}) == {-1,2}

Operations in matrix_multi<> take place in the folded subalgebra.

The common interface for framed_multi<> and matrix_multi<> is defined via the abstract base template class clifford_algebra<Scalar_T, Index_Set_T, Multivector_T>. This uses what I call pseudo-polymorphism. framed_multi<> and matrix_multi<> do not have a common base class, only a common interface defined via the template. The abstract base template class is used only to define the derived classes and is not used in the library. This keeps overhead to a minimum.

Non-member function code common between the two derived template classes is implemented via the use of template template functions.

For example, the functions sqrt, exp, log, cos, cosh, sin, sinh, tan, tanh, acos, acosh, asin, asinh, atan and atanh use algorithms based on matrix algorithms and Pade approximations, in such a way that the same algorithm is used for both framed_multi<> and matrix_multi<>.

The log and sqrt functions are based on Pade approximation and on matrix algorithms recently published by Cheng, Higham, Kenney and Laub.

The exp function is based on Pade approximation, and the other transcendental functions are based on exp and log. References for these are Abramowitz and Stegun, Gerald and Wheatley, Golub and van Loan. See details below.

Limitations

GluCat uses the C++ Standard Library and the Boost Library. Make sure you have these installed and working before attempting to use GluCat. Use the Boost download page to obtain the Boost Library.

Make sure that you have installed Boost with the same version of g++ as you will be using for GluCat.

The headers need a compiler which correctly handles template template classes.

GluCat 0.4.0 has not been tested with Boost versions before 1.33.0.

With g++ 4.0.0 and g++ 4.0.1, optimization with -O3 may corrupt iterators. Use -O3 -fno-strict-aliasing instead. See gcc bug 23599.

GluCat 0.4.0 has so far been tested using:

• Dual Intel(R) Pentium(R) 4, Kubuntu 7.10, gcc version 4.1.3, Boost 1.33.1
• AMD Athlon(tm) 64, SuSE Linux 10.1, with
• g++ 4.1.0, g++ 4.2.2, icpc 10.0 Build 20070426,
• Boost 1.33.0, 1.33.1, 1.34.0, 1.34.1.

To Do

There is plenty of testing, cleanup, packaging and documentation work yet to do.

• Improve the packaging of the example and test programs.
• Write a programmer's guide with descriptions of the API.
• Test and suggest improvements to John Fletcher's SWIG bindings to Ruby and Henk Jansen's SWIG bindings to Python.
• Simplify and improve the use of automake and autoconf in building GluCat.
• Investigate if GluCat could be used as the back-end to Christian Perwass' CLU.
• Test with Blitz++, deal.II, MTL, POOMA and uBLAS and adjust accordingly.
• Port to other architectures and compilers which support template template parameters.

Transcendental functions still need improvement.

• Fix sqrt() and log() so that they deal properly with negative real eigenvalues.
• Devise better algorithms and better implementations of existing algorithms for the transcendental functions. In particular, pay more attention to radius of convergence, condition number of matrices and poking out of the subalgebra.
• Try using the complex matrix representations. These use smaller matrices, which should be faster. This would mean changes to basis_element(), folded_dim() and others. Possibly deliver this as an additional template class rather than a change to matrix_multi<>.
• Try defining Boost concepts and more numeric traits so that GluCat can eventually become a Boost library.

The interface and library could be further improved in a number of different ways.

• Try refactoring the relationship between matrix_multi, framed_multi and clifford_algebra to allow more flexibility with template parameters. Possibly use enable_if and SFINAE to do this.
• Try adding a Matrix_Tag template parameter to framed_multi and matrix_multi, to determine if matrix_t is compressed, dense, etc.
• Try removing the template parameters LO and HI from framed_multi and matrix_multi, and using DEFAULT_LO and DEFAULT_HI where these are needed in framed_multi.h, matrx_multi.h, etc.
• Try making Tune_P into a mixin policy class template parameter.
• Add operator| for geometric operations: return rhs * lhs * involute(inv(rhs))
• Add convenience constructors to index_set<>: index_set<>(int, int), index_set<>(int, int, int), ... etc.
• Provide better support and testing for framed_multi< complex<> > etc.
• Try replacing multiplication by +/-1 within inner products by addition and subtraction.
• Try creating a class, vector_multi<>, which uses std::vector rather than std::map. This should be faster than framed_multi<>, if tuned properly, for full multivectors. For sparse multivectors, it may be slower.

Presentations

Publications

Authors

Paul C. Leopardi

Project supervisor

Dr. William McLean, School of Mathematics, UNSW.

Prototype

Arvind Raja

Arvind Raja's prototype was modelled on Clical, by Pertti Lounesto, R. Mikkola, V. Vierros.

A. Alexandrescu, N. Josuttis, A. Koenig, Scott Meyers, B. Moo, Bjarne Stroustrup, T. Veldhuizen:
C++ usage.
Joerg Arndt:
Bit wizardry.
C++, debugging.
Christian Perwass:
Operators, user interface.

Definitions

Milton Abramowicz and Irene A. Stegun:
Transcendental functions.
Mark Ashdown, G. P. Wene:
Negative integer indices.
Chris Doran, David Hestenes, Anthony Lasenby, Pertti Lounesto, Ian R. Porteous, Garret Sobczyk:
Clifford algebra sum, difference, inner product, outer product, geometric product, involutions and anti-involutions, norms.
Leo Dorst:
Left contraction.
Leo Dorst, Arvind Raja:
Index sets.
David Hestenes, Garret Sobczyk, G. P. Wene:
Countably infinite dimensional Clifford algebra.
Pertti Lounesto, Ian R. Porteous:
Matrix basis.

Algorithms

Milton Abramowicz and Irene A. Stegun, Gene H. Golub, Charles F. van Loan, C.F. Gerald, P.O. Wheatley:
Transcendental functions.
Joerg Arndt:
Bit functions on index sets.
Sheung Hun Cheng, Nicholas J. Higham, Charles S. Kenney, Alan J. Laub:
Product form of Denman-Beavers square root iteration, incomplete square root cascade logarithm.
Leo Dorst:
Left contraction.
P. Fleckenstein:
framed_multi<> geometric product
Gene H. Golub, Charles F. van Loan, Nicholas J. Higham:
Matrix division, iterative refinement.
Gene H. Golub, Charles F. van Loan:
Matrix exponential.
Gene H. Golub, Charles F. van Loan, C.F. Gerald, P.O. Wheatley, Doron Zeilberger:
Pertti Lounesto, Ian R. Porteous:
Matrix basis
.
A. Lumsdaine, J. Siek, MTL project:
Matrix algorithms and linear algebra.
Arvind Raja:
framed_multi<> sum, difference, inner product, outer product, geometric product, involutions and anti-involutions, norms.
Joerg Walter, uBLAS project:
LU factor and solve, matrix algorithms and linear algebra.
Jan Cnops:
Recursive expressions for matrix representations of Clifford algebras.
Michael Clausen, Ulrich Baum, David K. Maslen, Daniel N. Rockmore:
Generalized FFTs for finite groups.

Joerg Walter.

Configuration

Stephan Kulow, Walter Tasin, KDevelop team:
Configuration code used in config.in.in

Testing

Joerg Walter, John Fletcher, Henk Jansen, Johannes Brunen.

Related Projects

John Fletcher: BoostCliffordDiscussion

Workarounds

Carlos O'Ryan: workaround used in generator.h for: "... only defines a private destructor and has no friends"