Axiom/FriCAS - Grassmann/Clifford SPAD Code

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

• The full revised spad code is listed here.
• I have put the test results on this page.
• I have put user guide on this page.
• Requirements for further improvements on this page.

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
element	grade=m	grade = n-m

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

What I did was this:

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

• Added functions above.
• 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

• Fix the issues described at the end of this page.
• 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
          grade1 := gradeTerm(op1type) -- grade of first operand
          grade2 := gradeTerm(op2type) -- grade of second operand
          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)_
              +  _/_\(gradeInvolutionTerm(1$K,uType),lcProdTerm(1$K,op1type,1$K,vType)) -- gradeinvolution(u) /\ (x _| v)
            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))_
              +  _/_\(rcProdTerm(1$K,vType,1$K,op2type),gradeInvolutionTerm(1$K,uType)) -- (v |_ x)/\ gradeinvolution(u)
            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
          grade1 := gradeTerm(op1type) -- grade of first operand
          grade2 := gradeTerm(op2type) -- grade of second operand
          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
          grade1 := gradeTerm(op1type) -- grade of first operand
          grade2 := gradeTerm(op2type) -- grade of second operand
          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
• sign := (-1)^((gradeTerm(b1)-1)/2)

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()

Grade Involution

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

Code

        -- implements gradeInvolution for a single term by using:
        -- x = ((-1)^grade(x))*x
        gradeInvolutionTerm(mult: K,type1:SINT): % ==
          resul:% := New; resul.type1:=mult
          g:NNI := gradeTerm(type1)
          sign:Integer := (-1$Integer)^g
          resul := sign*resul
          resul

        -- implements gradeInvolution for a multivector by involution
        -- of each term seperately using:
        -- x = ((-1)^grade(x))*x
        gradeInvolution(x:%): % ==
            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

Grade

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

Code

        -- return the grade of the term, for example
        -- 1 is grade 0
        -- e1 is grade 1
        -- e1^e2 is grade 2 and so on
        gradeTerm(b:SINT): NNI ==
          grade:NNI := 0
          for i in 0..n-1 repeat
            if bit?(b,i) then grade := grade + 1
          grade

        -- return the max grade of multivector, for example
        -- 1 is grade 0
        -- e1 is grade 1
        -- e1^e2 is grade 2 and so on
        grade(x:%): NNI ==
          gr:NNI := 0         
          for ix in 0..dim-1 repeat
            if x.ix ~= 0 then
              gr := max(gr,gradeTerm(ix::SINT))
          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

• Grassman Algebra Spad Code
• CliffordPackage
• Clifford Algebra examples
• Clifford Algebra
• Clifford code

EuclideanSpace

home

This is one of a series of pages describing code that I have written for FriCAS. Some of this code has already been included with the FriCAS library and so can be used as soon as FriCAS is installed (as described here). Other code on these pages is not yet included with FriCAS and so would have to be compiled from within FriCAS. For an overview of all the code I have written, on various topics, see page here.

SPAD Language

This implementation of Grassmann/Clifford Algebra is written in the SPAD (a contraction of Scratchpad) language which is the user language of the Axiom/FriCAS computer algebra system.

For more information about the SPAD language see this page. If you would like to look at a simpler example of the language first see this page.

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

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.

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.

Axiom/Fricas Information on this Site

SPAD Language

Other Options - Idris