# Axiom/FriCAS - Grassmann/Clifford SPAD Code

I have modified the existing CliffordAlgebra domain in Axiom/FriCAS.

• The full revised spad code is listed here.

Dr Bertfried Fauser helped me with the theory.

In the longer term these algebras may be derived from the implementation that Bertfried is doing on Hopf algebras, as a (non Hopf) Hopf algebra deformation of the exterior algebra.

## Reason For Changes

When using this algebra for geometry and physics we want to be able to mix the Clifford product with the Grassmann (exterior or wedge) product and inner products in the same equation..

For example we may want to create a geometric object using meet and join then transform using conjugation.

In physics we need to combine these products by analogy with the way that vector algebra combines scalar, dot and cross products. When modeling solid bodies (isometries) we want to model either 'projective space' or 'conformal space'.

We also want these to work with a basis from non-diagonal quadratic form (such as bases which square to +ve and -ve and then rotated).

## Why So Many Product Types?

The exterior and inner products can be introduced in different ways, one way is to look at the exterior product as the product which generates the structure:

<e1, e2…en | ei/\ei=0, ei/\ej= -ej/\ei>

and the interior product defines the metric structure.

Or we can look at exterior and inner products as duals:

concept dual concept
product exterior inner
geometric interpretation

intersection of spaces

meet

union of spaces

join

Or we can look at the exterior and inner products as components of the Clifford product.

When we are dealing with pure vectors in orthogonal bases these concepts all coincide, however when we move to higher grades or use metrics with a non-diagonal quadratic form, the inner product types start to diverge and we end up with the regressive inner product, the contraction inner products, and so on.

It would be useful to have a complete table of all the products and their properties and uses here. Also it would be good to have a table for the various reverses, inverses and duals.

## Changes so Far

I have added the following functions (detailed notes on each of these at the end of this page):

 /\ _/_\: (%, %) -> % ++ Implement exterior Grassmann product operator ++ need to check precedence when used as an infix operator \/ _\_/: (%, %) -> % ++ Implement regressive inner,meet product operator ++ need to check precedence when used as an infix operator lc lc: (%, %) -> % ++ left contraction inner product rc rc: (%, %) -> % ++ right contraction inner product ~ _~: (%) -> % ++ reverse, complement,canonical dual basis gradeInvolution gradeInvolution: (%) -> % ++ implements gradeInvolution for a multivector by involution ++ of each term seperately using: ++ x = ((-1)^grade(x))*x eFromBinaryMap(NNI) eFromBinaryMap: NNI -> % ++ eFromBinaryMap(n) sets the appropriate unit elements, for example: ++ eFromBinaryMap(0) = 1 (scalar) ++ eFromBinaryMap(1) = e1 ++ eFromBinaryMap(2) = e2 ++ eFromBinaryMap(3) = e1*e2 ePseudoscalar() ePseudoscalar: () -> % ++ unit pseudoscalar grade(%) grade: (%) -> NNI ++ return the max grade of multivector, for example ++ 1 is grade 0 ++ e1 is grade 1 ++ e1^e2 is grade 2 and so on toTable() toTable: ((%, %) -> %) -> Matrix % ++ displays multiplication table for binary operation which ++ is represented as a function with two parameters. ++ row number represents first operand in binary order ++ column number represents second operand in binary order ++ could have returned type 'List List %' but matrix displays better toTable: ((%) -> %) -> Matrix % ++ displays table of unary function such as inverse, reverse, ++ complement, or dual basis ++ could have returned type 'List List %' but matrix displays better examples: toTable(/\) -- exterior table toTable(\/) -- regression table toTable(*) -- clifford table toTable(lc) -- left contraction table toTable(rc) -- right contraction table toTable(~) -- dual table

1. I edited the source file: fricas-1.0.8/src/algebra/clifford.spad.pamphlet
2. I extracted the spad code from )abbrev domain CLIF CliffordAlgebra to just before the @ symbol (Note: this does not include the QuadraticForm code which is not affected)
3. To avoid any name clashes I renamed all occurrences of CliffordAlgebra to be GrassmannAlgebra and CLIF to be GRAS
4. I saved this file as axiom/Grassmann.spad.
5. I added the functions above.

## Done So Far

• Modified the Clifford multiplication '*' so that it takes account of the non-diagonal terms in the bilinear form matrix.
• Changed constructor which previously required a QuadraticForm as a parameter. It now takes a SquareMatrix as input to represent a more general bilinear form. In other words there is no longer a requirement for the bilinear form to be symmetrical about the diagonal.
• Removed inheritance from VectorSpace, This domain still inherits from module but it does not yet support non-commutative elements.

## Further Work

• A grassmann basis for the clifford algebra (see Theoretical Questions above)
• A way to convert to and from multivectors and matrices (Pauli matrices)
• Support for equation solver for multivectors using Clifford, or possibly even mixed, products.
• Support for multivectors whose elements are non-commutative algebras.
• Allow multiple types of Clifford Product in the same multivector. Perhaps we could have an option of implementing the Clifford Product as a function which takes the two operands and also an alternative quadratic form.
• Move to implementation planned by Bertfried which is GrassmannSuperHopfAlgebraWithBasis
• Allow symbolic values for:
• indexes (e(i) instead of e(1)) which would be propagated to higher grade terms.
• bilinear forms with symbolic entries

I have put notes about implementing these things and more on the requirements page here.

## Exterior Product

Notation: /\

Returns 0 if there are any common terms, otherwise bitwise-or the bases and calculate the sign by putting both bases into a word and apply rule ei*ej = -ej*ei to get bases into numerical order.

Possible future extension:

• wedgeB:(%,%) -> R  (the extended bilinear form B : VxV -> R to B^ : /\Vx/\V ->R)

Code:

```        -- Implement exterior grassmann product on individual term in a
-- multivector and add it to the multivector.
exteriorProdTerm(op1mult: K, op1type:SINT, op2mult: K, op2type:SINT): % ==
resul:% := New
And(op1type,op2type) ~= 0 => resul -- if common terms return without adding
c := op1mult * op2mult
bz := Or(op1type,op2type)-- combine terms from both operands
for i in 0..n-1 | bit?(op1type,i) repeat
-- Apply rule  ei*ej = -ej*ei for i ~= j
k := 0
for j in i+1..n-1 | bit?(op1type, j) repeat k := k+1
for j in 0..i-1   | bit?(bz, j) repeat k := k+1
if odd? k then c := -c
resul.bz := resul.bz + c
resul

-- Implement exterior grassmann product
_/_\(x:%,y:%): % ==
z := New
for ix in 0..dim-1 repeat
if x.ix ~= 0 then for iy in 0..dim-1 repeat
if y.iy ~= 0 then z := z + exteriorProdTerm(x.ix,ix::SINT,y.iy,iy::SINT)
z```

## Regressive Inner Product

Notation: \/

This is implemented as the reversion of the exterior product:

~(x \/ y) = ~y /\ ~x

I have not implemented reversion yet (see below) so this is temporary code.

Temporary Code:

```        -- Implement regressive inner,meet product on individual term in a
-- multivector and add it to the multivector.
regressiveProdTerm(op1mult: K, op1type:SINT, op2mult: K, op2type:SINT): % ==
op1 := New; op1.op1type := 1\$K
op2 := New; op2.op2type := 1\$K
res := _/_\(op2*ePseudoscalar(),op1*ePseudoscalar()) *ePseudoscalar()
res := (op1mult * op2mult)*res
res

-- Implement regressive inner,meet product operator
_\_/(x:%,y:%): % ==
z := New
for ix in 0..dim-1 repeat
if x.ix ~= 0 then for iy in 0..dim-1 repeat
if y.iy ~= 0 then z := z + regressiveProdTerm(x.ix,ix::SINT,y.iy,iy::SINT)
z```

## Clifford Product

```        -- Implement Clifford multiplication on individual term in a
-- multivector.
cliffordProdTerm(op1mult: K, op1type:SINT, op2mult: K, op2type:SINT): % ==
resul:% := New
if grade1 = 0 then -- 1st operand is scalar so return scalar product
resul.op2type := resul.op2type + op1mult*op2mult
return resul
if grade2 = 0 then -- 2nd operand is scalar so return scalar product
resul.op1type := resul.op1type + op1mult*op2mult
return resul
if grade1 = 1 and grade2 = 1 then -- vect * vect = bilinear + x/\y
resul.(0\$NNI) := resul.(0\$NNI) + op1mult * op2mult * bilinear(op1type,op2type) -- add scalar term
resul := resul + exteriorProdTerm(op1mult,op1type,op2mult,op2type) -- add exterior term
return resul
if grade1 = 1 then -- 1st operand is vector so apply:
-- x*(u/\v) = x /\ u /\ v + x _| (u/\v)
-- = x/\u/\v + (x_|u)/\v + gradeinvolution(u) /\ (x _| v)
uType:SINT := leftMostBase(op2type) -- highest ^factor
vType:SINT := xor(op2type,uType) -- remaining ^factors
xt:% := New; xt.op1type := 1\$K
ut:% := New; ut.uType := 1\$K
vt:% := New; vt.vType := 1\$K
resul := _/_\(xt,exteriorProdTerm(1\$K,uType,1\$K,vType))_
+  _/_\(lcProdTerm(1\$K,op1type,1\$K,uType),vt)_
resul := (op1mult*op2mult)*resul -- apply 'scalar' multipliers
return resul
if grade2 = 1 then -- 2nd operand is vector so apply:
-- (v/\u)*x = v /\ u /\ x + rc(v/\u,x)
-- = v/\u/\x + v/\rc(u,x) +  rc(u, x)/\ gradeinvolution(v)
uType:SINT := rightMostBase(op1type) -- lowest ^factor
vType:SINT := xor(op1type,uType) -- remaining ^factors
xt:% := New; xt.op2type := 1\$K
ut:% := New; ut.uType := 1\$K
vt:% := New; vt.vType := 1\$K
resul := _/_\(exteriorProdTerm(1\$K,vType,1\$K,uType),xt)_
+  _/_\(vt,rcProdTerm(1\$K,uType,1\$K,op2type))_
resul := (op1mult*op2mult)*resul -- apply 'scalar' multipliers
return resul
-- if none of the above is true then apply:
-- (u /\ x) * v = u * (x _| v + x /\v) - (u |_ x) * v
xType:SINT := rightMostBase(op1type) -- highest ^factor
uType:SINT := xor(op1type,xType) -- remaining ^factors
ut:% := New; ut.uType := 1\$K
vt:% := New; vt.op2type := 1\$K
-- factor1 := x * v = x _| v + x /\v
-- factor1:% := cliffordProdTerm(1\$K,xType,1\$K,op2type)
factor1:% := lcProdTerm(1\$K,xType,1\$K,op2type)_
+ exteriorProdTerm(1\$K,xType,1\$K,op2type)
-- factor2 := u |_ x
factor2:% := rcProdTerm(1\$K,uType,1\$K,xType)
resul := ut * factor1 - factor2 * vt
if debug then
s1 := concat[toStringTerm(op1mult,op1type),"*",toStringTerm(op2mult,op2type)]
s2 := concat["=",toString(ut),"*",toString(factor1)]
s3 := concat["-",toString(factor2),"*",toString(vt)]
s4 := concat["=",toString(resul)]
sayTeX\$Lisp concat ["cliffordProdTerm: ",s1,s2,s3,s4]
resul := (op1mult*op2mult)*resul -- apply 'scalar' multipliers
resul```

## Left Contraction Inner Product

Notation: lc and rc (can't use _| because _ is escape character)

Code:

```        -- Implement left contraction on individual term in a
-- multivector and add it to the multivector.
lcProdTerm(op1mult: K, op1type:SINT, op2mult: K, op2type:SINT): % ==
resul:% := New
if grade1 = 0 then -- 1st operand is scalar so return scalar product
resul.op2type := resul.op2type + op1mult*op2mult
return resul
grade2 = 0 => resul -- 2nd operand is scalar so return 0
if grade1 = 1 and grade2 = 1 then -- vect _| vect = bilinear
resul.(0\$NNI) := resul.(0\$NNI) + op1mult * op2mult * bilinear(op1type,op2type) -- add scalar term
return resul
if grade1 = 1 then -- 1st operand is vector so apply: x_|(u^v)=(x_|u)^v - u^(x_|v)
uType:SINT := leftMostBase(op2type) -- highest ^factor
vType:SINT := xor(op2type,uType) -- remaining ^factors
inner2:% := New; inner2.vType := 1\$K
left:% := _/_\(lcProdTerm(op1mult, op1type,op2mult, uType),inner2)
inner4:% := New; inner4.uType := 1\$K
resul := resul + left + _/_\(inner4,lcProdTerm(-op1mult, op1type,op2mult, vType))
return resul
-- if none of the above is true then apply:
-- (u/\v) _| w = u _| (v _| w)
uType:SINT := leftMostBase(op1type) -- highest ^factor
vType:SINT := xor(op1type,uType) -- remaining ^factors
inner2:% := New
inner2.uType := 1\$K
resul := resul+ lc(inner2,lcProdTerm(op1mult,vType,op2mult,op2type))
resul

-- Implement left contraction inner product
lc(x:%,y:%): % ==
z := New
for ix in 0..dim-1 repeat
if x.ix ~= 0 then for iy in 0..dim-1 repeat
if y.iy ~= 0 then z := z + lcProdTerm(x.ix,ix::SINT,y.iy,iy::SINT)
z
```

## Right Contraction Inner Product

Notation: lc and rc (can't use _| because _ is escape character)

Code:

```        -- Implement right contraction on individual term in a
-- multivector and add it to the multivector.
rcProdTerm(op1mult: K, op1type:SINT, op2mult: K, op2type:SINT): % ==
resul:% := New
if grade2 = 0 then -- 1st operand is scalar so return scalar product
resul.op1type := resul.op1type + op1mult*op2mult
return resul
grade1 = 0 => resul -- 2nd operand is scalar so return 0
if grade1 = 1 and grade2 = 1 then -- vect |_ vect = bilinear
resul.(0\$NNI) := resul.(0\$NNI) + op1mult * op2mult * bilinear(op1type,op2type) -- add scalar term
return resul
if grade2 = 1 then -- 1st operand is vector so apply: (v^u)|_x=v^(u|_x) - (v|_x)^u
uType:SINT := rightMostBase(op1type) -- lowest ^factor
vType:SINT := xor(op1type,uType) -- remaining ^factors
inner2:% := New; inner2.vType := 1\$K
right:% := _/_\(inner2,rcProdTerm(op1mult, uType,op2mult, op2type))
inner4:% := New; inner4.uType := 1\$K
resul := resul + _/_\(rcProdTerm(op1mult, vType,-op2mult, op2type),inner4) + right
return resul
-- if none of the above is true then apply:
-- w |_ (v/\u) = (w |_ v) |_ u
uType:SINT := rightMostBase(op2type) -- lowest ^factor
vType:SINT := xor(op2type,uType) -- remaining ^factors
inner2:% := New
inner2.uType := 1\$K
resul := resul+ rc(rcProdTerm(op1mult,op1type,op2mult,vType),inner2)
resul

-- Implement right contraction inner product
rc(x:%,y:%): % ==
z := New
for ix in 0..dim-1 repeat
if x.ix ~= 0 then for iy in 0..dim-1 repeat
if y.iy ~= 0 then z := z + rcProdTerm(x.ix,ix::SINT,y.iy,iy::SINT)
z```

## Reverse

Notation: ~

The existing CliffordAlgebra also has inverse function:

recip(x: %): Union(%, "failed") -- Clifford inverse where possible

Currently the code only multiplies by the pseudoscalar, a lot more work to do on this!

If the basis is 'diagonal' then the reverse of an individual term could be calculated from:

• type := Xor(b1::SINT,(dim-1)::SINT) -- bitwise-not of binary map

Since this needs to cope with the 'non-diagonal' case I need to generalise the following examples. Bertfried has explained that this needs to split into 2 types:

Grassmann reversion Clifford reversion
~(e12) = e21 = -e12
~(e(12)) = ~(e(1)*e(2))
= ~(e12 +B(e1,e2))
= -e12+B(e1,e2)
= -e(1)*e(2) + 2B(e1,e2)
~(e(12)) = e(21) = e(2)*e(1) = e21+B(e2,e1)= -e12+B(e2,e1)
= -e(12)+B(e1,e2)+B(e2,e1)
now if B=B^t then we get
= -e(12)+2B(e1,e2)
if B=g+A with g=g^t and A=-A^t we get
= -e(12) +2g(e1,e2)~(e12)
= ~(e(12) - B(e1,e2))
= -e(12) +2B(e1,e2)-B(e1,e2)   (sym B)
= -e(12) + B(e1,e2)

Where:

• e12:=e1/\e2
• e(12):=e(1)*e(2) etc
• ~=reversion wrt the chosen basis

Temporary Code

```        -- Implement reverse, complement,canonical dual basis
_~(x:%): % ==
x * ePseudoscalar()
```

For each term seperately, change sign (- <-> +) if odd grade

Code

```        -- implements gradeInvolution for a single term by using:
resul:% := New; resul.type1:=mult
sign:Integer := (-1\$Integer)^g
resul := sign*resul
resul

-- implements gradeInvolution for a multivector by involution
-- of each term seperately using:
z := New
for ix in 0..dim-1 repeat
if x.ix ~= 0 then z := z + gradeInvolutionTerm(x.ix,ix::SINT)
z

```

## Constants

In addition to the e(1) unit vectors I have added constants for higher order bases.

• eFromBinaryMap(NNI)
• ePseudoscalar()

Note: CliffordAlgebra already has constants for vector bases and scalars:

• e(NNI)
• 1

Code

```        -- eFromBinaryMap(n) sets the appropriate unit elements, for example,
-- eFromBinaryMap(0) = 1  (scalar)
-- eFromBinaryMap(1) = e1
-- eFromBinaryMap(2) = e2
-- eFromBinaryMap(3) = e1*e2
eFromBinaryMap b ==
b >= dim => error "Too big"
z := New; z.b := 1; z

-- unit pseudoscalar
ePseudoscalar() ==
z := New
i := (dim-1)::NNI
z.i := 1
z
```

Returns the maximum grade of all the non-zero multivector terms.

Code

```        -- return the grade of the term, for example
-- e1^e2 is grade 2 and so on
for i in 0..n-1 repeat

-- return the max grade of multivector, for example
-- e1^e2 is grade 2 and so on
gr:NNI := 0
for ix in 0..dim-1 repeat
if x.ix ~= 0 then
gr
```

## Tables

Although these products are computed each time they are needed and not kept in a table (for the reasons explained above) I still find it useful to have a function which (when we are working in low dimensions) I can use to output the table as a matrix.

Code

```        -- displays multiplication table for binary operation which
-- is represented as a function with two parameters.
-- row number represents first operand in binary order
-- column number represents second operand in binary order
-- could have returned type 'List List %' but matrix displays better
toTable(fn:(%, %) -> %) ==
l: List % := [eFromBinaryMap(i) for i in 0..dim-1]
matrix [[fn(k,j) for j in l] for k in l]

-- displays table of unary function such as inverse, reverse, complement,
-- or dual basis
-- could have returned type 'List List %' but matrix displays better
toTable(fn:(%) -> %) ==
l: List % := [eFromBinaryMap(i) for i in 0..dim-1]
matrix [[j for j in l],[fn(k) for k in l]]
```

## Discussion

See these threads on FriCAS forum      Geometric Algebra for Physicists - This is intended for physicists so it soon gets onto relativity, spacetime, electrodynamcs, quantum theory, etc. However the introduction to Geometric Algebra and classical mechanics is useful.      Geometric Algebra for Computer Science: An Object-oriented Approach to Geometry. This book stresses the Geometry in Geometric Algebra, although it is still very mathematically orientated. Programmers using this book will need to have a lot of mathematical knowledge. Its good to have a Geometric Algebra book aimed at computer scientists rather than physicists. There is more information about this book here.      Clifford Algebras: Applications to Mathematics, Physics and Engineering (Progress in Mathematical Physics) (Hardcover) also see arXiv:math-ph/0212032 Clifford and Grassmann Hopf algebras via the BIGEBRA package for Maple Rafal Ablamowicz, B. Fauser: Comp. Physics Comm. 170, 2005:115--130 available from the arXiv