[GiNaC-list] Geometric Algebra

[Vladimir Kisil]

Alan Bromborsky (AB) wrote:
  AB>  To do this I must do a grade decomposition of multivectors. 

  I attach a sample program which defines a function to this end:

lst clifford_multi_decomp(const ex & e, const ex & c)

  For Clifford algebra in N dimensional space it returns a list of 2^n
coefficients in the standard Clifford algebra basis, i.e. element of the
list with the binary number 101 corresponds to e~2e~0, 110 -- to e~2e~1,
010 -- to e~1, etc. Please note, that the ordering of the units in the
poducts like e~2e~0 is determined by GiNaC in canonicalize_clifford()
procedure.

  Once clifford_multi_decomp() will be sufficiently tested I can prepare a
patch for its inclusion to the GiNaC core library.

      Best,
      	   Vladimir
--
Vladimir V. Kisil     email: kisilv at maths.leeds.ac.uk
--                      www: http://maths.leeds.ac.uk/~kisilv/

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

#define _ex2 numeric(2)

/* An auxiliary function used by simplify_indexed() and expand_dummy_sum() 
 * It returns an exvector of factors from the supplied product */
static void product_to_exvector(const ex & e, exvector & v, bool & non_commutative)
{
	// Remember whether the product was commutative or noncommutative
	// (because we chop it into factors and need to reassemble later)
	non_commutative = is_exactly_a<ncmul>(e);

	// Collect factors in an exvector, store squares twice
	v.reserve(e.nops() * 2);

	if (is_exactly_a<power>(e)) {
		// We only get called for simple squares, split a^2 -> a*a
		GINAC_ASSERT(e.op(1).is_equal(_ex2));
		v.push_back(e.op(0));
		v.push_back(e.op(0));
	} else {
		for (size_t i=0; i<e.nops(); i++) {
			ex f = e.op(i);
			if (is_exactly_a<power>(f) && f.op(1).is_equal(_ex2)) {
				v.push_back(f.op(0));
				v.push_back(f.op(0));
			} else if (is_exactly_a<ncmul>(f)) {
				// Noncommutative factor found, split it as well
				non_commutative = true; // everything becomes noncommutative, ncmul will sort out the commutative factors later
				for (size_t j=0; j<f.nops(); j++)
					v.push_back(f.op(j));
			} else
				v.push_back(f);
		}
	}
}

void product_to_blade(const ex & f, const ex & c, lst & res) {
	int  num = 0;
	// Collect factors in an exvector
	exvector v, s;
	unsigned char rl = ex_to<clifford>(c).get_representation_label();
	
	// Remember whether the product was commutative or noncommutative
	// (because we chop it into factors and need to reassemble later)
	bool non_commutative;
	product_to_exvector(f, v, non_commutative);

	exvector::const_iterator cit = v.begin(), citend = v.end();
	while (cit != citend) {
		if (is_a<clifford>(*cit) && is_a<cliffordunit>(cit->op(0)) && (ex_to<clifford>(*cit).get_representation_label() == rl) 
			&& ex_to<clifford>(*cit).same_metric(c)) {
			if (ex_to<varidx>(cit->op(1)).is_numeric())
				num = num + (int)exp2(ex_to<numeric>(ex_to<varidx>(cit->op(1)).get_value()).to_int());
			else
				throw(std::invalid_argument("product_toblade(): cannot decompose into blades with symbolic indices"));
		} else if (is_a<ncmul>(*cit))
			product_to_blade(*cit, c, res);
		else
			s.push_back(*cit);
		cit++;
	}
	res.let_op(num) = res.op(num) + remove_dirac_ONE(non_commutative ? ex(ncmul(s)) : ex(mul(s)), rl);
}


lst clifford_multi_decomp(const ex & e, const ex & c) {
	if(!ex_to<varidx>(c.op(1)).is_dim_numeric())
		throw(std::invalid_argument("clifford_multi_decom (): cannot decompose into blades with symbolic dimension"));

	int D = ex_to<numeric>(ex_to<varidx>(c.op(1)).get_dim()).to_int();
	lst res;
	for(size_t i=0; i < exp2(D); i++)
		res.append(0);

	ex e_exp = canonicalize_clifford(expand_dummy_sum(e, true));

	if (is_a<mul>(e_exp) || is_a<ncmul>(e_exp))
		product_to_blade(e_exp, c, res);
	else if (is_a<add>(e_exp))
		for (size_t i=0; i < e_exp.nops(); i++)
			product_to_blade(e_exp.op(i), c, res);
	else
		throw(std::invalid_argument("clifford_multi_decom (): cannot decompose into blades this object"));

	return res;
}

int main() {
	realsymbol a("a"), s("s");
	symbol x("x"), y("y");
	varidx mu(symbol("mu"), 3),  nu(symbol("nu"), 3);
	ex M = diag_matrix(lst(-1,-1,s));
	ex c = clifford_unit(mu, M, 1);

	ex e; 
	e = 2*clifford_unit(varidx(1, 3), M, 1)+x*clifford_unit(varidx(1, 3), M)*clifford_unit(varidx(2, 3, 1), M);
	cout << clifford_multi_decomp(e, c) << endl  << endl;
	e = 2*a*dirac_gamma(varidx(2,3))*dirac_gamma(varidx(1,3))*
		clifford_unit(varidx(2,3), M, 1)*clifford_unit(varidx(1,3), M, 1)
		+x*clifford_unit(mu, M, 1)*clifford_unit(mu.toggle_variance(), M, 1)
		+clifford_unit(varidx(0,3), M, 1)*clifford_unit(varidx(1,3), M)
		+y*clifford_unit(mu, M, 1)*indexed(matrix(1,3,lst(1,-1,1)),mu.toggle_variance())
		+5*clifford_unit(varidx(0,3), M, 1)*clifford_unit(varidx(1,3), M, 1)*clifford_unit(varidx(2,3), M, 1);
	cout << clifford_multi_decomp(e, c) << endl  << endl;

}