Axiom/FriCAS - Grassmann/Clifford SPAD Code

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

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

Further Work

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:

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:

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:

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.

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

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


metadata block
see also:
Correspondence about this page

Book Shop - Further reading.

Where I can, I have put links to Amazon for books that are relevant to the subject, click on the appropriate country flag to get more details of the book or to buy it from them.

flag flag flag flag flag flag 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.
flag flag flag flag flag flag 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.
flag flag flag flag flag flagClifford 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

 

This site may have errors. Don't use for critical systems.

Copyright (c) 1998-2017 Martin John Baker - All rights reserved - privacy policy.