Generic library of universal Clifford algebra templates (original) (raw)

GluCat: Generic library of universal Clifford algebra templates

$Real representation matrix for neutral Clifford algebra R_{2,2}$	SourceForge.net project page Introduction Presentations Publications Authors Home page

SourceForge.net project page GluCat project pagehosted at:

About GluCat

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 asBlitz++,deal.II,Matrix Template Library andPOOMA. 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 _framed_multivector, 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 inGA 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 functionseven, frame, inv, involute, reverse, pure,quad, scalar etc.

The aim is to get the notation as close as possible to that used byfloat 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...HIas 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. 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 onstd::map<index_set<LO,HI>,Scalar_T> or (optionally) on __gnu_cxx::hash_map< const index_set<LO,HI>,Scalar_T, hash <LO,HI> >orstd::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 theBoost 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 classgenerator_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 classclifford_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 atanhuse 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 onmatrix 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

Paul Leopardi,"A generalized FFT for Clifford algebras", Bulletin of the Belgian Mathematical Society - Simon Stevin, Volume 11, Number 5, 2005, pp. 663-688.
Preprint:Applied Maths Report AMR04/17, March, 2004.

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.

Advice

A. Alexandrescu, N. Josuttis, A. Koenig, Scott Meyers, B. Moo, Bjarne Stroustrup, T. Veldhuizen:

C++ usage.

Joerg Arndt:

Bit wizardry.

Duraid Madina, Russell Standish:

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:

Pade approximation.

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.

uBLAS Interface

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.

John Fletcher:BoostCliffordDiscussion

Workarounds

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

References

Jörg Arndt,Algorithms for programmers,Chapter 1, Bit Wizardry, as of 2007-02-28.
Milton Abramowicz and Irene A. Stegun, Handbook of mathematical functions, Dover 1972, first published 1965.
Sheung Hun Cheng, Nicholas J. Higham, Charles S. Kenney, Alan J. Laub,Approximating The Logarithm Of A Matrix To Specified Accuracy, 1999.
Michael Clausen, Ulrich Baum, Fast Fourier transforms, Wissenschaftsverlag, Manheim, 1993.
Leo Dorst, "Honing geometric algebra for its use in the computer sciences", pp. 127-152 of G. Sommer, (ed.)Geometric Computing with Clifford Algebras, Springer 2001.
C.F. Gerald and P.O. Wheatley, Applied Numerical Analysis, 6th Edition, Addison-Wesley, 1999.
Gene H. Golub and Charles F. van Loan, Matrix Computations, 3rd ed., Johns Hopkins UP, 1996.
Nicholas J. Higham,Accuracy and Stability of Numerical Algorithms, SIAM, 1996.
Nicholas J. Higham, "The Scaling and Squaring Method for the Matrix Exponential Revisited", SIAM Journal on Matrix Analysis and Applications, Vol. 26, Issue 4 (2005), pp. 1179-1193.
Pertti Lounesto,"Clifford algebra calculations with a microcomputer", pp. 39-55 of A. Micali, R. Boudet, J. Helmstetter, (eds),Clifford algebras and their applications in mathematical physics : proceedings of second workshop held at Montpellier, France, 1989, Dordrecht ; Boston : Kluwer Academic Publishers, 1992.
Pertti Lounesto,Clifford Algebras and Spinors, Cambridge University Press, 2001.
Ian R. Porteous,Clifford Algebras and the Classical GroupsCambridge University Press, 1995.
Arvind Raja,"Object-oriented implementations of Clifford algebras in C++: a prototype", inClifford Algebras with Numeric and Symbolic Computations, Rafal Ablamowicz, P. Lounesto, Josep M. Parra (Editors), 1996.
G. P. Wene,"The Idempotent stucture of an infinite dimensional Clifford algebra", pp. 161-164 ofClifford algebras and their applications in mathematical physics, 1992, as above.
Doron Zeilberger,"PADE" (Maple code), 2002.

Home page

Home page of Paul Leopardi at ANU

Page Generated: Wed 7 November 2007