
Boost : 
From: Jaap Suter (boost_at_[hidden])
Date: 20040610 00:26:12
Hi,
[crossposted to oonumerics and boostdevelopers]
About 20 months back, I asked if there was any interest in a geometric or
Clifford algebra library written in C++. Despite some positive replies, I
got the
impression that there were not a lot of people in the Boost community who
saw any practical uses for such library.
Back then I had a rudimentary proofofconcept for a Clifford algebra
library that used metaprogramming techniques to gain significant
performance benefits over other existing libraries (Gaigen, nKlein, CLU,
etcetera).
Over the past few months I have put a lot more work into this library. The
bulk of the work consisted of adding more unittests, proving many of the
theorems in some of the papers by David Hestenes et al. I also added
expression templates to improve performance.
To give an idea of what the library is capable of:
#include <boost/math/clifford/clifford.hpp>
#include <boost/mpl/list_c.hpp>
#include <cassert>
using namespace boost::math::clifford;
using namespace boost::mpl;
typedef algebra<4, 1> conformal_model;
typedef blade<float, 0, conformal_model> scalar;
typedef blade<float, 1, conformal_model> vector;
typedef blade<float, 5, conformal_model> pseudo_scalar;
typedef blade<float, 2, conformal_model> bi_vector;
// same as:
// typedef element<float, list_c<int, 1>, conformal_model> bi_vector;
typedef list_c<int, 0, 2> spinor_grades;
typedef element<float, spinor_grades, conformal_model> spinor;
typedef list_c<int, 0, 1, 2, 3, 4, 5> multi_vector_grades;
typedef element<float, multi_vector_grades, conformal_model> multi_vector;
int main()
{
assert(sizeof(scalar) == sizeof(float));
assert(sizeof(vector) == 5 * sizeof(float));
assert(sizeof(bi_vector) == 10 * sizeof(float));
assert(sizeof(pseudo_scalar) == sizeof(float));
assert(sizeof(spinor) == 11 * sizeof(float));
assert(sizeof(multi_vector) == 32 * sizeof(float));
vector e_plus;
e_plus[3] = 1.0;
vector e_min;
e_min[4] = 1.0;
e_star = e_min  e_plus;
e_star *= std::sqrt(2.0);
e_star /= 2.0f;
e [noalias] = e_plus * std::sqrt(2.0)  e_star;
bi_vector E = outer_product(e, e_star);
// geometric product between bi_vector and vector
multi_vector v = E * e_star;
// geometric product between two vectors, generating a spinor
vector a, b;
spinor s = a * b;
// doing a spinor rotation
vector c;
c [noalias] = reverse(s) * a * s;
}
Each product always does the bare minimum of additions, multiplications and
subtractions to calculate the result on the lefthandside of an expression.
Thus, if you do the full geometric product of two multi_vectors but assign
it to a scalar, it will only do those calculations that are involved in
grade 0.
Also, unlike other existing libraries, no dynamic allocations are being done
ever.
Currently, the following operations are supported, between any two random
elements of an algebra (blades and composite multi_vectors):
geometric_product = operator *
inner_product
outer_product
fat_dot_product
left_contraction
right_contraction
scalar_product
of course, all the normal operations (scalar multiplication, addition of
elements, etcetera) are supported, as well as reversion, involution and
conjugates.
Adding more operators is trivial. They can either be implemented in terms of
existing operators, or by implementing an appropriate gradepredicate. For
example, this is what the implementation of the outer_product looks like:
struct outer_product_predicate
{
template <class ResultGrade, class LhsGrade, class RhsGrade>
struct apply
{
typedef typename mpl::if_
<
typename mpl::equal_to
<
ResultGrade,
typename mpl::plus<LhsGrade,RhsGrade>::type
>::type,
mpl::true_,
mpl::false_
>::type type;
};
};
BOOST_MATH_CLIFFORD_PRODUCT_IMPL(
outer_product,
impl::outer_product_predicate
)
which for somebody that knows the definition of the outer_product is easy to
understand.
I have unittests for 2, 3 4 and 5 dimensional algebras, but more is
possible if your compiler supports it.
Currently I have only tested the library on MSVC 7.1. I tried 7.0 but
couldn't get it to work (excessive use of MPL). GCC and Intel compiled it at
one point, but because I only use MSVC I dropped support. Of course, if
anybody is interested I'd be trying other compilers again. You will need a
compiler that does aggresive inlining though. On MSVC I am using the
force_inline feature to get the release performance I need.
Anyway, let me know if anybody is interested, and I'll put together a
package to upload somewhere.
Cheers,
Jaap
Boost list run by bdawes at acm.org, gregod at cs.rpi.edu, cpdaniel at pacbell.net, john at johnmaddock.co.uk