clifford.spad.pamphlet

\documentclass{article}
\usepackage{axiom}
\usepackage{amssymb}
% \input{diagrams.tex}
\begin{document}
\title{\$SPAD/src/algebra clifford.spad}
\author{Martin J Baker, Stephen M. Watt}
\maketitle
\begin{abstract}
\end{abstract}
\eject
\tableofcontents
\eject
\section{domain QFORM QuadraticForm}
<<domain QFORM QuadraticForm>>=
)abbrev domain QFORM QuadraticForm
++ Author: Stephen M. Watt
++ Date Created: August 1988
++ Date Last Updated: May 17, 1991
++ Basic Operations: quadraticForm, elt
++ Related Domains: Matrix, SquareMatrix
++ Also See:
++ AMS Classifications:
++ Keywords: quadratic form
++ Examples:
++ References:
++
++ Description:
++   This domain provides modest support for quadratic forms.
QuadraticForm(n, K): T == Impl where
    n: PositiveInteger
    K: Field
    SM ==> SquareMatrix
    V  ==> DirectProduct

    T ==> AbelianGroup with
        quadraticForm: SM(n, K) -> %
            ++ quadraticForm(m) creates a quadratic form from a symmetric,
            ++ square matrix m.
        matrix: % -> SM(n, K)
            ++ matrix(qf) creates a square matrix from the quadratic form qf.
        elt: (%, V(n, K)) -> K
            ++ elt(qf,v) evaluates the quadratic form qf on the vector v,
            ++ producing a scalar.
        if SM(n, K) has ConvertibleTo InputForm then ConvertibleTo InputForm


    Impl ==> SM(n,K) add
        Rep := SM(n,K)

        quadraticForm m ==
            not symmetric? m =>
                error "quadraticForm requires a symmetric matrix"
            m::%
        matrix q == q pretend SM(n,K)
        elt(q,v) == dot(v, (matrix q * v))

        if SM(n, K) has ConvertibleTo InputForm then
            convert(q:%):InputForm ==
               -- without package call we get infinite recursion
               mif := (convert(matrix q)$(SM(n, K)))@InputForm
               qf := convert(convert('quadraticForm)@SExpression)@InputForm
               convert([qf, mif])$InputForm


@
\section{domain CLIF CliffordAlgebra\cite{7,12}}
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).
\subsection{Why So Many Product Types?}
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:
\begin{list}{}
topic 	concept 	dual concept
product 	exterior 	inner
geometric interpretation 	intersection of spaces	meet
	union of spaces	join
element 	grade=m 	grade = n-m
\end{list}

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. 

\subsection{instantiate an instance of CliffordAlgebra}
We first need to instantiate an instance of CliffordAlgebra, the constructor
takes 3 parameters:
\begin{list}{}
parameter 	type 	meaning
1st 	PositiveInteger 	dimension of the underling algebra (module over)
2nd 	Type 	type over which this is defined
3rd 	SquareMatrix of dimension defined by 1st parameter 	bilinear form
\end{list}

In the following example we will define a 3 dimensional Euclidean space so
the bilinear form is the identity matrix:

(1) -> Eu := CliffordAlgebra(3,Fraction(Integer),[[1,0,0],[0,1,0],[0,0,1]])                

   (1)  CliffordAlgebra(3,Fraction(Integer),[[1,0,0],[0,1,0],[0,0,1]])
                                                            Type: Domain
\subsection{Entering Multivectors}
Entering Multivectors

We can then create some elements in the algebra (I will call the general
element in this algebra a 'multivector' to be consistent with most of the
literature). First we can create some basis vectors using e(PI):

(2) -> e1:Eu := e(1)                                                         

   (2)  e
         1
      Type: CliffordAlgebra(3,Fraction(Integer),[[1,0,0],[0,1,0],[0,0,1]])
(3) -> e2:Eu := e(2)                                                         

   (3)  e
         2
      Type: CliffordAlgebra(3,Fraction(Integer),[[1,0,0],[0,1,0],[0,0,1]])

There are many options for creating higher grade basis,

    * by multiplying basis vectors,
    * or by using ee(List PI) to specify the higher order term in one go
      with a list of its indices (exterior product of indicies)
    * or from a binary map of the indicies
    * or using the monomial function (Clifford product of indicies)
    * or using the pseudoscalar.

(4) -> e12:Eu := e(1) /\ e(2)

   (4)  e e
         1 2
      Type: CliffordAlgebra(3,Fraction(Integer),[[1,0,0],[0,1,0],[0,0,1]])
(5) -> e12:Eu := ee[1,2]

   (5)  e e
         1 2
      Type: CliffordAlgebra(3,Fraction(Integer),[[1,0,0],[0,1,0],[0,0,1]])
(6) -> e12:Eu := eFromBinaryMap(3)

   (6)  e e
         1 2
      Type: CliffordAlgebra(3,Fraction(Integer),[[1,0,0],[0,1,0],[0,0,1]])
(7) -> e12:Eu := monomial(1,[1,2])

   (7)  e e
         1 2
      Type: CliffordAlgebra(3,Fraction(Integer),[[1,0,0],[0,1,0],[0,0,1]])
(8) -> p:Eu := ePseudoscalar()

   (8)  e e e
         1 2 3
      Type: CliffordAlgebra(3,Fraction(Integer),[[1,0,0],[0,1,0],[0,0,1]])

We can create vectors from a linear sum of its bases and similarly we can
create a bivector from a sum of bases with grade 2:

(9) -> v:Eu := 2*e(1) + 3*e(2) + 4*e(3)

   (9)  2e  + 3e  + 4e
          1     2     3
        Type: CliffordAlgebra(3,Fraction(Integer),[[1,0,0],[0,1,0],[0,0,1]])
(10) -> m:Eu := 2*ee[1,2] + 3*ee[1,3] + 4*ee[2,3]

   (10)  2e e  + 3e e  + 4e e
          1 2     1 3     2 3
        Type: CliffordAlgebra(3,Fraction(Integer),[[1,0,0],[0,1,0],[0,0,1]])

If we want a more complicated multivector, then it may be easier to specify
it in one, by using the multivector function. Note that the indices are
specified in binary order,
that is multivector[scalar,e1,e2,e12,e3,e13,e23,e123...] so we only need
specify the multipliers for these:

(11) -> m:Eu := multivector[0,0,0,2,0,3,4,0]

   (11)  2e e  + 3e e  + 4e e
          1 2     1 3     2 3
        Type: CliffordAlgebra(3,Fraction(Integer),[[1,0,0],[0,1,0],[0,0,1]])

\subsection{Algebraic Operations}
Algebraic Operations

We can then do some operations such as:
\begin{list}{}
operator 	meaning
/\ 	exterior
\/ 	regression
* 	Clifford
lc 	left contraction
rc 	right contraction
\end{list}

examples:

(12) -> e1 /\ e2                                                              

   (12)  e e
         1 2
      Type: CliffordAlgebra(3,Fraction(Integer),[[1,0,0],[0,1,0],[0,0,1]])
(13) -> e1 /\ e1                                                              
                                                                                                    
   (13)  0                                                                                           
      Type: CliffordAlgebra(3,Fraction(Integer),[[1,0,0],[0,1,0],[0,0,1]])                       
(14) -> e1 \/ e2                                                                                     
                                                                                                    
   (14)  0                                                                                           
      Type: CliffordAlgebra(3,Fraction(Integer),[[1,0,0],[0,1,0],[0,0,1]])                       
(15) -> e1 \/ e1                                                                                     
                                                                                                    
   (15)  0                                                                                           
      Type: CliffordAlgebra(3,Fraction(Integer),[[1,0,0],[0,1,0],[0,0,1]])                       
(16) -> lc(e1,e2)                                                                                    
                                                                                                    
   (16)  0                                                                                           
      Type: CliffordAlgebra(3,Fraction(Integer),[[1,0,0],[0,1,0],[0,0,1]])                       
(17) -> lc(e1,e1)                                                                                    
                                                                                                    
   (17)  1
      Type: CliffordAlgebra(3,Fraction(Integer),[[1,0,0],[0,1,0],[0,0,1]])
(18) -> rc(e1,e2)

   (18)  0
      Type: CliffordAlgebra(3,Fraction(Integer),[[1,0,0],[0,1,0],[0,0,1]])
(19) -> rc(e1,e1)

   (19)  1
      Type: CliffordAlgebra(3,Fraction(Integer),[[1,0,0],[0,1,0],[0,0,1]])

\subsection{Functions}
Functions

Here are some functions available for the algebra :
\begin{list}{}
coefficient(%,List NNI) 	extracts the scalar coefficient of a given basis which is specified by a list of indices.
~(%) 	Clifford dual
gradeInvolution(%) 	invert (multiply by -1) vector and modify other grades accordingly
reverse(%) 	reverse order of basis vectors in each blade
conj(%) 	equivalent to applying both gradeInvolution and reverse
grade(%) 	returns an integer representing the maximum grade of the non-zero blades in this multivector. 0 for scalar, 1 for vector and so on.
\end{list}

examples:

(20) -> c:Eu := coefficient(p,[1,2])

   (20)  0
      Type: CliffordAlgebra(3,Fraction(Integer),[[1,0,0],[0,1,0],[0,0,1]])
(21) -> c:Eu := coefficient(p,[1,2,3])

   (21)  1
      Type: CliffordAlgebra(3,Fraction(Integer),[[1,0,0],[0,1,0],[0,0,1]])
(22) -> ~e12

   (22)  - 1
      Type: CliffordAlgebra(3,Fraction(Integer),[[1,0,0],[0,1,0],[0,0,1]])
(23) -> gradeInvolution(e12)

   (23)  e e
         1 2
      Type: CliffordAlgebra(3,Fraction(Integer),[[1,0,0],[0,1,0],[0,0,1]])
(24) -> gradeInvolution(p)

   (24)  - e e e
           1 2 3
      Type: CliffordAlgebra(3,Fraction(Integer),[[1,0,0],[0,1,0],[0,0,1]])
(25) -> grade(p)

   (25)  3
                                                      Type: PositiveInteger
\subsection{Tables}
Tables

We can produce tables (Cayley table) of the built-in operations, unary and
binary functions by using 'toTable'. The result is a matrix showing all
combinations of the basis for each operand:

(26) -> toTable(lc)$Eu

         +1  e   e   e e   e   e e   e e   e e e +
         |    1   2   1 2   3   1 3   2 3   1 2 3|
         |                                       |
         |0  1   0    e    0    e     0     e e  |
         |             2         3           2 3 |
         |                                       |
         |0  0   1   - e   0    0     e    - e e |
         |              1              3      1 3|
         |                                       |
         |0  0   0   - 1   0    0     0     - e  |
   (26)  |                                     3 |
         |                                       |
         |0  0   0    0    1   - e   - e    e e  |
         |                        1     2    1 2 |
         |                                       |
         |0  0   0    0    0   - 1    0      e   |
         |                                    2  |
         |                                       |
         |0  0   0    0    0    0    - 1    - e  |
         |                                     1 |
         |                                       |
         +0  0   0    0    0    0     0     - 1  +
Type: Matrix(CliffordAlgebra(3,Fraction(Integer),[[1,0,0],[0,1,0],[0,0,1]]))
(27) -> toTable(*)$Eu

         +  1       e        e       e e      e       e e      e e    e e e +
         |           1        2       1 2      3       1 3      2 3    1 2 3|
         |                                                                  |
         |  e       1       e e       e      e e       e      e e e    e e  |
         |   1               1 2       2      1 3       3      1 2 3    2 3 |
         |                                                                  |
         |  e     - e e      1       - e     e e    - e e e     e     - e e |
         |   2       1 2                1     2 3      1 2 3     3       1 3|
         |                                                                  |
         | e e     - e       e       - 1    e e e    - e e     e e     - e  |
         |  1 2       2       1              1 2 3      2 3     1 3       3 |
   (27)  |                                                                  |
         |  e     - e e    - e e    e e e     1       - e      - e     e e  |
         |   3       1 3      2 3    1 2 3               1        2     1 2 |
         |                                                                  |
         | e e     - e    - e e e    e e      e       - 1     - e e     e   |
         |  1 3       3      1 2 3    2 3      1                 1 2     2  |
         |                                                                  |
         | e e    e e e     - e     - e e     e       e e      - 1     - e  |
         |  2 3    1 2 3       3       1 3     2       1 2                1 |
         |                                                                  |
         |e e e    e e     - e e     - e     e e       e       - e     - 1  |
         + 1 2 3    2 3       1 3       3     1 2       2         1         +
Type: Matrix(CliffordAlgebra(3,Fraction(Integer),[[1,0,0],[0,1,0],[0,0,1]]))

\subsection{Debugging}
Debugging

There are some commands to help with debugging:

setMode("orthogonal",false)

The "orthogonal" flag is set to true if the non-diagonal terms in the matrix
are zero. This allows a more optimised Clifford multiplication for this
case. The above setMode command would turn that optimisation off forcing all
multiplications to be done by the more general algorithm. This allows the
output of the two algorithms to be compared.

setMode("debug",true)

This "debug" flag causes the internal working of the algorithms to be
displayed during calculations.

<<domain CLIF CliffordAlgebra>>=
)abbrev domain CLIF CliffordAlgebra
++ Author: Martin Baker (adapted from a version by Stephen M. Watt)
++ Date Created: August 1988
++ Date Last Updated: Feb 2010
++ Basic Operations: wholeRadix, fractRadix, wholeRagits, fractRagits
++ Related Domains: QuadraticForm, Quaternion, Complex
++ Also See:
++ AMS Classifications:
++ Keywords: clifford algebra, grassmann algebra, spin algebra
++ Examples:
++ References:
++
++ Description:
++  CliffordAlgebra(n, K, bLin) defines a module of dimension \spad{2^n}
++  over K, given a bilinear form bLin on \spad{K^n}.
++
++  Examples of Clifford Algebras are: gaussians, quaternions, exterior
++  algebras and spin algebras.

CliffordAlgebra(n, K, bLin): T == Impl where
    n: PositiveInteger
    K: Field
    bLin: SquareMatrix(n, K)

    PI ==> PositiveInteger
    NNI==> NonNegativeInteger
    SINT==> SingleInteger

    T ==> Join(Ring, Algebra(K)) with
        e: PI -> %
          ++ e(n) produces phi(e_i) where e_i is i-th basis
          ++ vector in K^n and phi is canonical embedding of
          ++ K^n into Clifford algebra.
        ee: List PI -> %
          ++ to allow entries like: ee[1,2]
        multivector: List K -> %
          ++ to allow entries like: 1+2*e1+3*e2+4*e1e2 = multivector[1,2,3,4]
        coerce: % -> List(K)
          ++ coerce to multivector
        coerce: List(K) -> %
          ++ coerce from multivector (same as multivector function above)
        eFromBinaryMap: NNI -> %
          ++ eFromBinaryMap(n) sets the appropriate Grassmann basis, for
          ++ example:
          ++ eFromBinaryMap(0) = 1  (scalar)
          ++ eFromBinaryMap(1) = e1
          ++ eFromBinaryMap(2) = e2
          ++ eFromBinaryMap(3) = e1/\e2
        ePseudoscalar: () -> %
          ++ unit pseudoscalar
        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
        monomial: (K, List PI) -> %
          ++ monomial(c,[i1,i2,...,iN]) produces the value given by
          ++ \spad{c*e(i1)*e(i2)*...*e(iN)}.
        coefficient:  (%, List PI) -> K
          ++ coefficient(x,[i1,i2,...,iN])  extracts the coefficient of
          ++ \spad{e(i1)*e(i2)*...*e(iN)} in x.
        recip: % -> Union(%, "failed")
          ++ recip(x) computes the multiplicative inverse of x or "failed"
          ++ if x is not invertible.
        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
        _/_\: (%, %) -> %
          ++ Implement exterior grassmann product operator
          ++ need to check precidence when used as an infix operator
        _\_/: (%, %) -> %
          ++ Implement regressive inner,meet product operator
          ++ need to check precidence when used as an infix operator
        lc: (%, %) -> %
          ++ left contraction inner product
        rc: (%, %) -> %
          ++ right contraction inner product
        _~: (%) -> %
          ++ reverse, complement,canonical dual basis
        gradeInvolution: (%) -> %
          ++ x = ((-1)^grade(x))*x
        reverse: (%) -> %
          ++ implements reverse for a single term by using:
          ++ grade: 0  1  2  3...
          ++ multi: 1  1 -1 -1...
        conj: (%) -> %
          ++ implements Clifford conjugate for a multivector by involution
          ++ and reverse of each term seperately using:
          ++ grade: 0  1  2  3...
          ++ multi: 1 -1 -1  1...
        setMode: (String,Boolean) -> Boolean
          ++ allows override of parameters such as orthogonal
          ++ used for debugging


    Impl ==> add
        Qeelist:Vector K :=  diagonal(bLin)::Vector K
          ++ contains the diagonal terms of bLin (what the base vectors
          ++ square to)
        orthogonal:Boolean :=  diagonal?(bLin)
          ++ true if bininear form is diagonal (non diagonal terms are zero)
        debug:Boolean := false
          ++ if set to true displays working when calculating results
        dim     :=  2^n
          ++ the total dimension of all the grades from scalar up to
          ++ pseudoscalar
        Rep     := PrimitiveArray K
          ++ array of Field which can hold the multivector values
        New     ==> new(dim, 0$K)$Rep

        x, y, z: %
        c: K
        m: Integer

        characteristic() == characteristic()$K
        dimension(): CardinalNumber == dim::CardinalNumber

        x = y ==
            for i in 0..dim-1 repeat
                if x.i ~= y.i then return false
            true

        x + y == (z := New; for i in 0..dim-1 repeat z.i := x.i + y.i; z)
        x - y == (z := New; for i in 0..dim-1 repeat z.i := x.i - y.i; z)
        - x   == (z := New; for i in 0..dim-1 repeat z.i := - x.i; z)
        m * x == (z := New; for i in 0..dim-1 repeat z.i := m*x.i; z)
        c * x == (z := New; for i in 0..dim-1 repeat z.i := c*x.i; z)
        -- we want module over non-commutative
        -- but no available signatures for:
        --x * m == (z := New; for i in 0..dim-1 repeat z.i := x*m.i; z)
        --x * c == (z := New; for i in 0..dim-1 repeat z.i := x*c.i; z)

        0            == New
        1            == (z := New; z.0 := 1; z)
        coerce(m): % == (z := New; z.0 := m::K; z)
        coerce(c): % == (z := New; z.0 := c; z)

        -- output single term to string
        -- this produces a simple one line string unlike coerceMonom
        -- no formatting such as subscript is done
        -- this is intended mainly for debug purposes
        toStringTerm(c: K, b: SINT): String ==
            mult:String := (mathObject2String$Lisp c)@String
            b = 0 => mult
            c = 0 => "0"
            if c=1 then mult := "e"
            if c= -1 then mult := "-e"
            if c~=1 and c~= -1 then mult := concat(mult,"e")
            for i in 0..n-1 repeat
              if bit?(b,i) then mult := concat(mult,string(i+1)$String)
            mult

        -- output multivector to string
        -- this produces a simple one line string unlike coerce(x): Ex
        -- no formatting such as subscript is done
        -- this is intended mainly for debug purposes
        toString(m:%): String ==
            res:String := ""
            for im in 0..dim-1 repeat
              if m.im ~= 0 then
                s:String := (mathObject2String$Lisp m.im)@String
                neg:Boolean := false
                if #s > 0 then neg := (s.1 = char "-")
                if res ~= "" and not neg then
                  res := concat[res,"+",toStringTerm(m.im,im::SINT)]
                if res = "" or neg then
                  res := concat(res,toStringTerm(m.im,im::SINT))
            if res = "" then res:="0"
            res

        -- e(n) produces the appropriate unit element.
        --        iz   z
        -- e(1) = 1    e1
        -- e(2) = 2    e2
        -- e(3) = 4    e3
        e b ==
            b::NNI > n => error "No such basis element"
            iz := 2^((b-1)::NNI)
            z := New; z.iz := 1; z

        -- to allow entries like e[1,2]
        ee(l) == 
          lst:List % := [e i for i in l]
          reduce(_/_\, lst , 1)

        -- to allow entries like: 1+2*e1+3*e2+4*e1e2 = multivector[1,2,3,4]
        multivector(mv:List K) ==
          (z := New; for i in 0..dim-1 repeat z.i := mv(i+1); z) 

        -- convert multivector (list of multipliers of bases)
        -- to a \spadtype{CliffordAlgebra} representation
        coerce(mv:List(K)): % ==
          (z := New; for i in 0..dim-1 repeat z.i := mv(i+1); z)

        -- convert \spadtype{CliffordAlgebra} representation to a
        -- multivector (list of multipliers of bases).
        coerce(rp:%): List(K) == [rp.i for i in 0..dim-1]

        -- 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() ==
            p := New
            i := (dim-1)::NNI
            p.i := 1
            p

        -- 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]]

        -- 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 + e1 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

        -- 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 := if odd?(g) then -1 else 1
          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

        -- implements reverse for a single term by using:
        -- grade: 0  1  2  3...
        -- multi: 1  1 -1 -1...
        reverseTerm(mult: K,type1:SINT): % ==
          resul:% := New; resul.type1:=mult
          g := gradeTerm(type1)::SINT
          sign:Integer := if odd?(shift(g,-1)) then -1 else 1
          resul := sign*resul
          resul

        -- implements reversal for a multivector by reverse
        -- of each term seperately using:
        -- grade: 0  1  2  3...
        -- multi: 1  1 -1 -1...
        reverse(x:%): % ==
            z := New
            for ix in 0..dim-1 repeat
                if x.ix ~= 0 then z := z + reverseTerm(x.ix,ix::SINT)
            z

        -- implements Clifford conjugation for a single term by using:
        -- grade: 0  1  2  3...
        -- multi: 1 -1 -1  1...
        conjTerm(mult: K,type1:SINT): % ==
          resul:% := New; resul.type1:=mult
          g := gradeTerm(type1)::SINT
          sign:Integer := if odd?(shift(g+1,-1)) then -1 else 1
          resul := sign*resul
          resul

        -- implements Clifford conjugate for a multivector by involution
        -- and reverse of each term seperately using:
        -- grade: 0  1  2  3...
        -- multi: 1 -1 -1  1...
        conj(x:%): % ==
            z := New
            for ix in 0..dim-1 repeat
                if x.ix ~= 0 then z := z + conjTerm(x.ix,ix::SINT)
            z

        -- allows override of parameters such as orthogonal
        -- used for debugging
        setMode(s:String,val:Boolean): Boolean ==
          if s = "orthogonal" then (orthogonal := val; return true)
          if s = "debug" then (debug := val; return true)
          false

        -- return the base vector number from binary, for example
        -- binary for e1 is 001 which returns 1
        -- binary for e2 is 010 which  returns 2 and so on
        baseVect(b:SINT): NNI ==
          for i in 0..n-1 repeat
            if bit?(b,i) then return i+1
          0

        -- return a term from bilinear product
        bilinear(b1:SINT,b2:SINT): K ==
          bv1 := baseVect(b1)
          bv2 := baseVect(b2)
          if bv1 = 0$NNI then return 0
          if bv2 = 0$NNI then return 0
          bLin(bv1,bv2)

        -- for a given term, return the base furthest on the left, for
        -- example e1^e2^e3 would return e1
        leftMostBase(b:SINT): SINT ==
          mask:SINT := 1
          for i in 0..n-1 repeat
            if And(mask,b) ~= 0 then return mask
            mask := shift(mask,1)::SINT
          0

        -- for a given term, return the base furthest on the right, for
        -- example e1^e2^e3 would return e3
        rightMostBase(b:SINT): SINT ==
          mask:SINT := shift(1,(n-1)::SINT)$SINT
          for i in 0..n-1 repeat
            if And(mask,b) ~= 0 then return mask
            mask := shift(mask,-1)::SINT
          0

        -- Implement exterior grassmann product on individual term in a
        -- multivector.
        exteriorProdTerm(op1mult: K, op1type:SINT,_
                         op2mult: K, op2type:SINT): % ==
          resul:% := New
          -- if common terms return without adding
          And(op1type,op2type) ~= 0 => resul
          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 regressive inner,meet product on individual term in a
        -- 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 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

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

        -- Implement left contraction on individual term in a
        -- 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
            -- add scalar term
            resul.(0$NNI) := resul.(0$NNI) +_
                  op1mult * op2mult * bilinear(op1type,op2type)
            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 right contraction on individual term in a
        -- 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
            -- 2nd operand is scalar so return scalar product
            resul.op1type := resul.op1type + op1mult*op2mult
            return resul
          grade1 = 0 => resul -- 1st 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
            -- 2nd 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
            if debug then
              s1 := concat[toStringTerm(op1mult,op1type),"L",_
                                  toStringTerm(op2mult,op2type)]
              s2 := concat["= ",toStringTerm(op1mult, vType),"L",_
                                  toStringTerm(-op2mult, op2type)]
              s3 := concat["/\",toString(inner4),"+",toString(right)]
              s4 := concat("=",toString(resul))
              sayTeX$Lisp concat ["rcProdTerm: ",s1,s2,s3,s4]
            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)
          if debug then
            s1 := concat[toStringTerm(op1mult,op1type),"L",_
                         toStringTerm(op2mult,op2type)]
            s2 := concat["= (",toStringTerm(op1mult,op1type),"L",_
                         toStringTerm(op2mult,vType)]
            s3 := concat(") L",toString(inner2))
            s4 := concat("=",toString(resul))
            sayTeX$Lisp concat ["rcProdTerm: ",s1,s2,s3,s4]
          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

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

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

        -- Implement Clifford multiplication for orthogonal bases
        -- on individual term in a multivector.
        -- If the bilinear form is diagonal then this may be
        -- more efficient than cliffordProdTerm as it does not
        -- require recursion.
        -- The ei*ej products could instead be precomputed in
        -- a (2^n)^2 multiplication table although only 
        -- practical for low dimension size.
        cliffordDiagonalTerm(op1mult: K, op1type:SINT,_
                             op2mult: K, op2type:SINT): % ==
            c  := op1mult * op2mult
            bz := op2type
            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
                -- Apply rule  ei^2 = Q(ei)
                if bit?(bz,i) then
                    c := c * Qeelist.(i+1)
                    bz:= bz - (2^i)::SINT
                else
                    bz:= bz + (2^i)::SINT
            result := New; result.bz := c
            result

        -- Clifford product
        -- if orthogonal then this operation will use faster algorithm
        x * y ==
          z := New
          if orthogonal then
            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 + cliffordDiagonalTerm(x.ix,ix::SINT,y.iy,iy::SINT)
          if not orthogonal then
            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 + cliffordProdTerm(x.ix,ix::SINT,y.iy,iy::SINT)
          z

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

        -- used by monomial which is used by recip
        -- this only uses the diagonal terms in the bilinear form
        -- so may not be completely general
        canonMonom(c: K, lb: List PI): Record(coef: K, basel: NNI) ==
            -- 0. Check input
            for b in lb repeat b > n => error "No such basis element"
            -- 1. Apply identity ei*ej = -ej*ei, i ~= j.
            -- The Rep assumes n is small so bubble sort is ok.
            -- Using bubble sort keeps the exchange info obvious.
            wasordered   := false
            exchanges := 0
            while not wasordered repeat
                wasordered := true
                for i in 1..#lb-1 repeat
                    if lb.i > lb.(i+1) then
                        t := lb.i; lb.i := lb.(i+1); lb.(i+1) := t
                        exchanges := exchanges + 1
                        wasordered := false
            if odd? exchanges then c := -c
            -- 2. Prepare the basis element
            -- Apply identity ei*ei = Q(ei).
            bz := 0
            for b in lb repeat
                bn := (b-1)::NNI
                if bit?(bz, bn) then
                    c := c * Qeelist bn
                    bz:= ( bz - 2^bn )::NNI
                else
                    bz:= bz + 2^bn
            [c, bz::NNI]

        -- monomial(c,[i1,i2,...,iN]) produces the value given by
        -- \spad{c*e(i1)*e(i2)*...*e(iN)}.
        -- this uses canonMonom which only uses the diagonal terms in the
        -- bilinear form so may not be completely general
        -- (K, List PI) -> %
        monomial(c, lb) ==
            r := canonMonom(c, lb)
            z := New
            z r.basel := r.coef
            z

        -- coefficient(x,[i1,i2,...,iN])  extracts the coefficient of
        -- \spad{e(i1)*e(i2)*...*e(iN)} in x.
        -- (%, List PI) -> K
        coefficient(z, lb) ==
            r := canonMonom(1, lb)
            r.coef = 0 => error "Cannot take coef of 0"
            z r.basel/r.coef

        -- use Ex instead of OutputForm
        Ex ==> OutputForm

        -- output single term to output device
        -- vector bases have the form e1, e2 ...
        -- bivector bases have the form e12, e23 ...
        -- and so on
        coerceMonom(c: K, b: NNI): Ex ==
            b = 0 => c::Ex
            ml := [sub("e"::Ex, i::Ex) for i in 1..n | bit?(b,i-1)]
            be := reduce("*", ml)
            c = 1 => be
            c::Ex * be

        -- output multivector to output device
        -- as sum of terms, like,
        -- 1 + 2*e1 + 3*e2 + 4*e1e2
        coerce(x): Ex ==
            tl := [coerceMonom(x.i,i) for i in 0..dim-1 | x.i ~= 0]
            null tl => "0"::Ex
            reduce("+", tl)

        localPowerSets(j:NNI): List(List(PI)) ==
          l: List List PI := list []
          j = 0 => l
          Sm := localPowerSets((j-1)::NNI)
          Sn: List List PI := []
          for x in Sm repeat Sn := cons(cons(j pretend PI, x),Sn)
          append(Sn, Sm)

        powerSets(j:NNI):List List PI == map(reverse, localPowerSets j)

        Pn:List List PI := powerSets(n)

        -- Clifford multiplicative inverse
        -- not guaranteed to exist, if we cant
        -- find an inverse then we return the string "failed"
        -- this generates a matrix representation and then inverts matrix
        -- this uses monomial which uses canonMonom which only uses the
        -- diagonal terms in the bilinear form so may not be completely
        -- general 
        recip(x: %): Union(%, "failed") ==
          one:% := 1
          -- tmp:c := x*yC - 1$C
          rhsEqs : List K := []
          lhsEqs: List List K := []
          lhsEqi: List K
          for pi in Pn repeat
            rhsEqs := cons(coefficient(one, pi), rhsEqs)
            lhsEqi := []
            for pj in Pn repeat
                lhsEqi := cons(coefficient(x*monomial(1,pj),pi),lhsEqi)
            lhsEqs := cons(reverse(lhsEqi),lhsEqs)
          ans := particularSolution(matrix(lhsEqs),
            vector(rhsEqs))$LinearSystemMatrixPackage(K, Vector K, Vector K,_
                           Matrix K)
          ans case "failed" => "failed"
          ansP := parts(ans)
          ansC:% := 0
          for pj in Pn repeat
            cj:= first ansP
            ansP := rest ansP
            ansC := ansC + cj*monomial(1,pj)
          ansC

@
\section{domain PAULI PauliMatrix}
(1) -> )read axiom/testPauli.input
-- First we will create a Clifford algebra and its corresponding
-- Pauli matrix. The Clifford algebra is created from the bilinear form
-- and the Pauli matrix is created from basis matrices. At this stage,
-- we have to know the equivalents already, since I don't yet know an
-- algorithmn to convert between bilinear form and basis matrices.
)library CLIF

   CliffordAlgebra is now explicitly exposed in frame frame1
   CliffordAlgebra will be automatically loaded when needed from
      /home/martin/CLIF.NRLIB/CLIF
)library PAULI

   PauliMatrix is now explicitly exposed in frame frame1
   PauliMatrix will be automatically loaded when needed from
      /home/martin/PAULI.NRLIB/PAULI
K := Fraction(Integer)


   (1)  Fraction(Integer)
                                                                   Type: Type
CA := CliffordAlgebra(2,K,[[1,0],[0,1]])


   (2)  CliffordAlgebra(2,Fraction(Integer),[[1,0],[0,1]])
                                                                   Type: Type
PM := PauliMatrix(2,false,2,K,[matrix[[1,0],[0,-1]],matrix[[0,1],[1,0]]])


   (3)  PauliMatrix(2,false,2,Fraction(Integer),[MATRIX,MATRIX])
                                                                   Type: Type
--
-- we can get confidence that these two are equivalent by checking that
-- their multiplication tables are equivalent:
toTable(*)$CA


        + 1      e      e    e e +
        |         1      2    1 2|
        |                        |
        | e      1     e e    e  |
        |  1            1 2    2 |
   (4)  |                        |
        | e    - e e    1    - e |
        |  2      1 2           1|
        |                        |
        |e e    - e     e    - 1 |
        + 1 2      2     1       +
             Type: Matrix(CliffordAlgebra(2,Fraction(Integer),[[1,0],[0,1]]))
toTable()$PM


        + +1  0+    +1   0 +    +0  1+    + 0   1+ +
        | |    |    |      |    |    |    |      | |
        | +0  1+    +0  - 1+    +1  0+    +- 1  0+ |
        |                                          |
        |+1   0 +    +1  0+    + 0   1+    +0  1+  |
        ||      |    |    |    |      |    |    |  |
        |+0  - 1+    +0  1+    +- 1  0+    +1  0+  |
   (5)  |                                          |
        | +0  1+    +0  - 1+    +1  0+    +- 1  0+ |
        | |    |    |      |    |    |    |      | |
        | +1  0+    +1   0 +    +0  1+    + 0   1+ |
        |                                          |
        |+ 0   1+  + 0   - 1+  +1   0 +  +- 1   0 +|
        ||      |  |        |  |      |  |        ||
        ++- 1  0+  +- 1   0 +  +0  - 1+  + 0   - 1++
       Type: Matrix(PauliMatrix(2,false,2,Fraction(Integer),[MATRIX,MATRIX]))
--
-- we can coerse both ways between these domains but we have to coerce via
-- List. I can't write a direct coercion because, at compile time, the
-- Clifford algebra does not know the exact signature of Pauli matrix and
-- Pauli matrix does not know the exact signature of Clifford algebra.
-- Is there a way round this?
-- However we can see here that the round trip c1 -> p1 -> c2 gets us
-- back to where we started.
c1 := [1::K,2::K,3::K,4::K]::CA

   (6)  1 + 2e  + 3e  + 4e e
              1     2     1 2
                     Type: CliffordAlgebra(2,Fraction(Integer),[[1,0],[0,1]])
p1 := c1::List K::PM

        + 3    7 +
   (7)  |        |
        +- 1  - 1+
               Type: PauliMatrix(2,false,2,Fraction(Integer),[MATRIX,MATRIX])
p1 := c1::PM


   Cannot convert from type CliffordAlgebra(2,Fraction(Integer),[[1,0],
      [0,1]]) to PauliMatrix(2,false,2,Fraction(Integer),[MATRIX,MATRIX
      ]) for value
   1 + 2e  + 3e  + 4e e
         1     2     1 2

c2 := p1::List K::CA

   (8)  1 + 2e  + 3e  + 4e e
              1     2     1 2
                     Type: CliffordAlgebra(2,Fraction(Integer),[[1,0],[0,1]])
--
-- we can check that multipication is the same for non-basis values:
c3 := [5::K,6::K,7::K,8::K]::CA

   (9)  5 + 6e  + 7e  + 8e e
              1     2     1 2
                     Type: CliffordAlgebra(2,Fraction(Integer),[[1,0],[0,1]])
p3 := c3::List K::PM

         +11   15 +
   (10)  |        |
         +- 1  - 1+
               Type: PauliMatrix(2,false,2,Fraction(Integer),[MATRIX,MATRIX])
c1*c3

   (11)  6 + 20e  + 14e  + 24e e
                1      2      1 2
                     Type: CliffordAlgebra(2,Fraction(Integer),[[1,0],[0,1]])
p4 := p1*p3

         + 26    38 +
   (12)  |          |
         +- 10  - 14+
               Type: PauliMatrix(2,false,2,Fraction(Integer),[MATRIX,MATRIX])
c4 := p4::List K::CA

   (13)  6 + 20e  + 14e  + 24e e
                1      2      1 2
                     Type: CliffordAlgebra(2,Fraction(Integer),[[1,0],[0,1]])
--
-- we can also check that division is the same (at least for orthogonal basis)
p5 := [3/10::K,0::K,0::K,2/5::K]::PM

         + 3   2 +
         |--   - |
         |10   5 |
   (14)  |       |
         |  2   3|
         |- -  --|
         +  5  10+
               Type: PauliMatrix(2,false,2,Fraction(Integer),[MATRIX,MATRIX])
p6 := inverse(p5)$PM

         +6    8+
         |-  - -|
         |5    5|
   (15)  |      |
         |8   6 |
         |-   - |
         +5   5 +
    Type: Union(PauliMatrix(2,false,2,Fraction(Integer),[MATRIX,MATRIX]),...)
c5 := p5::List K::CA

          3   2
   (16)  -- + - e e
         10   5  1 2
                     Type: CliffordAlgebra(2,Fraction(Integer),[[1,0],[0,1]])
c6 := recip(c5)$CA

         6   8
   (17)  - - - e e
         5   5  1 2
          Type: Union(CliffordAlgebra(2,Fraction(Integer),[[1,0],[0,1]]),...)
c6 := p6::List K::CA

         6   8
   (18)  - - - e e
         5   5  1 2
                     Type: CliffordAlgebra(2,Fraction(Integer),[[1,0],[0,1]])
(19) ->

So far we have only considered an algebra based on a 2D space. In this case
the dimension of the multivector is 2^2 and the dimension of the Pauli matrix
is 2 by 2. Although the code is written for the 'n' dimensional case I don't
know if this is valid for 'n' dimensional space? Can we
always represent an 2^n dimensional multivector by a nxn dimensional Pauli
matrix? Usually 3D space is represented by a matrix of complex numbers, so
we have a multivector of dimension 2^3 represented by a 2x2 Pauli matrix
containing complex numbers. In what circumstances can we represent an 2^n
dimensional multivector by a (n-1)x(n-1) dimensional Pauli matrix?

How can I allow for both possibilities in the same domain? My first thought
was to specify the type and dimenstion of the multivector and matrix
seperately in the constructor. However, even if they are both of type Field,
we can't coerce or even 'pretend' between them. So then I thought of having a
boolean value cmplx in the constructor and then defining the types like this:

K : Field -- type of multivector
C ==> (if cmplx then Complex K else K) -- type of elements of matrix

but I can't use C as if it were were a type? I thought of using:

if cmplx then C : Complex K else C : K

but cmplx is only known at runtime so this won't compile. I even thought of
defining somthing like this:

C : Union(K,Complex K)

but I can't work how to use it?


<<domain PAULI PauliMatrix>>=
)abbrev domain PAULI PauliMatrix
++ Author: Martin J Baker
++ Date Created: Feb 2010
++ Date Last Updated: Feb 2010
++ Basic Operations:
++ Related Domains: SquareMatrix, CliffordAlgebra
++ Also See:
++ AMS Classifications:
++ Keywords: matrix, linear algebra, Pauli matrix, Clifford Algebra
++ Examples:
++ References:
++ Description:
++ The \spad{PauliMatrix} is a square matrix which can be represented by
++ a linear combination of basis matrices.
++ I'm not sure if this is the best name as I think this is a bit more
++ general than the matrices that arise in Pauli's treatment of spin
++ in quantum mechanics?
++ one application of this is to find matrix representation isomorphic
++ to a given Clifford Algebra
PauliMatrix(ndim,cmplx,mdim,K,Bases): Exports == Implementation where
  ndim : PositiveInteger
    ++ Dimension of Pauli matrix, that is a ndim by ndim SquareMatrix
  cmplx : Boolean
    ++ if true then matrix elements are Complex K if false then K
  mdim : PositiveInteger
    ++ Number of dimensions of underlying space. The dimension of the
    ++ multivector (or Clifford numbers if you prefer that terminology) will
    ++ therefore be mdim^2
    ++ To make more general why not define over Ring but need for
    ++ compatibility with Clifford Algebra stops me doing this.
  K : Field
    ++ Type of elements of multivector
    ++ Used Field for compatibility with Clifford Algebra
--  C ==> (if cmplx then Complex K else K)
--  if cmplx then C : Complex K else C : K
--  C ==> Union(K,Complex K)
  C ==> K
    ++ Type of elements of Pauli matrix, often complex numbers
    ++ To make more general why not define over Ring but use of Complex
    ++ (reqires at least CommutativeRing) stops me doing this.

  Bases: List(SquareMatrix(ndim,C))

  PI ==> PositiveInteger
  NNI==> NonNegativeInteger
  Row ==> DirectProduct(ndim,C)
  Col ==> DirectProduct(ndim,C)
  MATLIN ==> MatrixLinearAlgebraFunctions(C,Row,Col,%)

  Exports ==> Join(SquareMatrixCategory(ndim,C,Row,Col),_
                   CoercibleTo Matrix C) with
    transpose: % -> %
      ++ \spad{transpose(m)} returns the transpose of the matrix m.
    pauliMatrix: Matrix C -> %
      ++ \spad{pauliMatrix(m)} converts a matrix of type \spadtype{Matrix}
      ++ to a matrix of type \spadtype{pauliMatrix}.
    coerce: % -> Matrix C
      ++ \spad{coerce(m)} converts a matrix of type \spadtype{PauliMatrix}
      ++ to a matrix of type \spadtype{Matrix}.
    matrixToBigVector: SquareMatrix(ndim,C) -> Vector C
      ++ \spad{matrixToBigVector(m)} stacks the rows of the matrix m into
      ++ a big vector
    listColToMatrix: (List Vector C) -> SquareMatrix(ndim*ndim,C)
      ++ \spad{listColToMatrix(m)} combines a list of column vectors 
      ++ into a matrix
    coerce: % -> List(K)
    coerce: List(K) -> %
      ++ convert multivector (list of multipliers of bases)
      ++ to a \spadtype{PauliMatrix}.
      ++ I would have liked to coerce directly from PauliMatrix like:
      ++ coerce: CliffordAlgebra(...) -> %
      ++ but to do this we need to know the Type that we are converting from
      ++ at compile time but it is only known at runtime.
    eFromBinaryMap: NNI -> SquareMatrix(ndim,C)
      ++ eFromBinaryMap(n) sets the appropriate Grassmann basis, for example:
      ++ eFromBinaryMap(0) = 1  (scalar)
      ++ eFromBinaryMap(1) = e1
      ++ eFromBinaryMap(2) = e2
      ++ eFromBinaryMap(3) = e1/\e2
    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
    if C has CommutativeStar and C has unitsKnown then unitsKnown
      ++ the invertible matrices are simply the matrices whose determinants
      ++ are units in the Ring R.
    if C has ConvertibleTo InputForm then ConvertibleTo InputForm

  Implementation ==> Matrix C add

    Qelt2 ==> QAREF2O$Lisp
    Qsetelt2 ==> QSETAREF2O$Lisp

    qelt(m:%, i:Integer, j:Integer) == Qelt2(m, i, j, 1@Integer, 1@Integer)
    qsetelt_!(m:%,i:Integer,j:Integer,r:C):C ==
      Qsetelt2(m, i, j, r, 1@Integer, 1@Integer)

    scalarMatrix r == scalarMatrix(ndim,r)$Matrix(C) pretend %
    0    == scalarMatrix 0
    1    == scalarMatrix 1

    New:% ==> (new(ndim,ndim,0)$Matrix(C)) pretend %
      ++ macro to create a new representation

    -- \spad{matrixToBigVector(m)} stacks the rows of the matrix m into
    -- a big vector
    matrixToBigVector(m) ==
      bv :Vector(C) := new(ndim*ndim,0)
      for i in 0..(ndim-1) repeat
        for j in 0..(ndim-1) repeat
          bv.(i*ndim + j +1) := m(i+1,j+1)
      bv

    -- \spad{listColToMatrix(m)} combines a list of column vectors
    -- into a matrix
    listColToMatrix(lv) ==
      m1 := 0$SquareMatrix(ndim*ndim,C)
      for i in 1..(ndim*ndim) repeat
        setColumn!(m1 pretend Matrix C,i,lv.i)
      m1

    buildBaseMatrix:SquareMatrix(ndim*ndim,C) ==
      a:Vector(C) := new(ndim*ndim,0$C)$Vector(C)
      lv:List(Vector(C)) := [copy a for i in 1..(ndim*ndim)]
      for i in 0..((ndim*ndim)-1) repeat
        baseMat := eFromBinaryMap i
        v:Vector(C) := matrixToBigVector(baseMat)
        setelt(lv,i+1,v)
      listColToMatrix(lv)

    baseMatrix:SquareMatrix(ndim*ndim,C) := buildBaseMatrix()

    -- convert multivector (list of multipliers of bases)
    -- to a \spadtype{PauliMatrix}.
    -- I would have liked to coerce directly to PauliMatrix like:
    -- coerce: CliffordAlgebra(...) -> %
    -- but to do this we need to know the Type that we are converting from
    -- at compile time but it is only known at runtime.
    coerce(mv:List(K)): % ==
      -- first convert multvector from type K to type C
      mv2:List(C) := [(mv.i)::C for i in 1..(mdim*mdim)]
      bigVec:DirectProduct(ndim*ndim,C) := _
        baseMatrix * directProduct(vector(mv2))$DirectProduct(ndim*ndim,C)
      resu:% := New
      for i in 0..(ndim-1) repeat
        for j in 0..(ndim-1) repeat
          qsetelt_!(resu,i+1,j+1,(bigVec.(i*ndim + j +1))::C)
      resu

    -- convert \spadtype{PauliMatrix} to a
    -- multivector (list of multipliers of bases).
    -- I would have liked to coerce directly from PauliMatrix like:
    -- coerce: CliffordAlgebra(...) -> %
    -- but to do this we need to know the Type that we are converting from
    -- at compile time but it is only known at runtime.
    coerce(mv:%): List(K) ==
      lhsEqs := baseMatrix
      rhsEqs:Vector(C) := matrixToBigVector(mv pretend SquareMatrix(ndim,C))
      ans := particularSolution(lhsEqs,rhsEqs)$LinearSystemMatrixPackage(C,_
        Vector C, Vector C, SquareMatrix(ndim*ndim,C))
      ans case "failed" => 
        sayTeX$Lisp "failed"
        new(ndim*ndim,0)
      -- finally convert multvector from type C to type K
      ansList:List(C) := parts(ans)
      result:List(K) := [(ansList.i)::K for i in 1..(mdim*mdim)]
      result

    pauliMatrix x ==
      (nrows(x) ~= ndim) or (ncols(x) ~= ndim) =>
        error "pauliMatrix: matrix of bad dimensions"
      copy(x) pretend %

    coerce(x:%):Matrix(C) == copy(x pretend Matrix(C))
    row(x,i)    == directProduct row(x pretend Matrix(C),i)
    column(x,j) == directProduct column(x pretend Matrix(C),j)
    coerce(x:%):OutputForm == coerce(x pretend Matrix C)$Matrix(C)

    -- 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 >= mdim*mdim => error "Too big"
      baseMat := 1$SquareMatrix(ndim,C)
      for j in 1..(#Bases) repeat
        if bit?(b,j-1) then
          baseMat := baseMat * Bases.j
      baseMat

    -- 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() ==
      l: List % := [eFromBinaryMap(i) pretend % for i in 0..(mdim*mdim-1)]
      matrix([[k * j for j in l] for k in l])$Matrix(%)

    if C has CommutativeRing then
      determinant x == determinant(x pretend Matrix(C))
      minordet x    == minordet(x pretend Matrix(C))
      Pfaffian x    == Pfaffian(x pretend Matrix(C))

    if C has EuclideanDomain then
      rowEchelon x == rowEchelon(x pretend Matrix(C)) pretend %

    if C has IntegralDomain then
      rank x    == rank(x pretend Matrix(C))
      nullity x == nullity(x pretend Matrix(C))
      nullSpace x ==
        [directProduct c for c in nullSpace(x pretend Matrix(C))]

    if C has Field then
      dimension() == (m * n) :: CardinalNumber

      inverse x ==
        (u := inverse(x pretend Matrix(C))) case "failed" => "failed"
        (u :: Matrix(C)) pretend %

      x:% ^ n:Integer ==
        ((x pretend Matrix(C)) ^ n) pretend %

      recip x == inverse x

    if C has ConvertibleTo InputForm then
      convert(x:%):InputForm ==
         convert [convert('pauliMatrix)@InputForm,
                  convert(x::Matrix(C))]$List(InputForm)



@
\section{License}
<<license>>=
--Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
--All rights reserved.
--
--Redistribution and use in source and binary forms, with or without
--modification, are permitted provided that the following conditions are
--met:
--
--    - Redistributions of source code must retain the above copyright
--      notice, this list of conditions and the following disclaimer.
--
--    - Redistributions in binary form must reproduce the above copyright
--      notice, this list of conditions and the following disclaimer in
--      the documentation and/or other materials provided with the
--      distribution.
--
--    - Neither the name of The Numerical ALgorithms Group Ltd. nor the
--      names of its contributors may be used to endorse or promote products
--      derived from this software without specific prior written permission.
--
--THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
--IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
--TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
--PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
--OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
--EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
--PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
--PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
--LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
--NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
--SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
@
<<*>>=
<<license>>

<<domain QFORM QuadraticForm>>
<<domain CLIF CliffordAlgebra>>
@
\eject
\begin{thebibliography}{99}
\bibitem{1} Lounesto, P.
"Clifford algebras and spinors",
2nd edition, Cambridge University Press (2001)
\bibitem{2} Porteous, I.,
"Clifford algebras and the classical groups",
Cambridge University Press (1995)
Van Nostrand Reinhold, (1969)
\bibitem{3} Bergdolt, G.
"Orthonormal basis sets in Clifford algebras",
in \cite{16} (1996)
\bibitem{4} Dorst, Leo,
"Honing geometric algebra for its use in the computer sciences",
pp127-152 from \cite{15} (2001)
\bibitem{5} Braden, H.W.,
"N-dimensional spinors: Their properties in terms of finite groups",
American Institute of Physics,
J. Math. Phys. 26(4), April 1985
\bibitem{6} Lam, T.Y. and Smith, Tara L.,
"On the Clifford-Littlewood-Eckmann groups: a new look at periodicity mod 8",
Rocky Mountains Journal of Mathematics, vol 19, no. 3, (Summer 1989)
\bibitem{7} Leopardi, Paul "Quick Introduction to Clifford Algebras"\\
{\bf http://web.maths.unsw.edu.au/~leopardi/clifford-2003-06-05.pdf}
\bibitem{8} Cartan, Elie and Study, Eduard
"Nombres Complexes",
Encyclopaedia Sciences Math\'ematique, \'edition fran\c caise, 15, (1908),
d'apr\`es l'article allemand de Eduard Study, pp329-468. Reproduced as
pp107-246 of \cite{17}
\bibitem{9} Hestenes, David and Sobczyck, Garret
"Clifford algebra to geometric calculus: a unified language for
mathematics and physics", D. Reidel, (1984)
\bibitem{10} Wene, G.P.,
"The Idempotent structure of an infinite dimensional Clifford algebra",
pp161-164 of \cite{13} (1995)
\bibitem{11} Ashdown, M.
"GA Package for Maple V",\\
http://www.mrao.cam.ac.uk/~clifford/software/GA/GAhelp5.html
\bibitem{12} Doran, Chris and Lasenby, Anthony,
"Geometric Algebra for Physicists"
Cambridge University Press (2003) ISBN 0-521-48022-1
\bibitem{13} Micali, A., Boudet, R., Helmstetter, J. (eds),
"Clifford algebras and their applications in mathematical physics:
proceedings of second workshop held at Montpellier, France, 1989",
Kluwer Academic Publishers (1992)
\bibitem{14} Porteous, I.,
"Topological geometry"
Van Nostrand Reinhold, (1969)
\bibitem{15} Sommer, G. (editor),
"Geometric Computing with Clifford Algebras",
Springer, (2001)
\bibitem{16} Ablamowicz, R., Lounesto, P., Parra, J.M. (eds)
"Clifford algebras with numeric and symbolic computations",
Birkh\"auser (1996)
\bibitem{17} Cartan, Elie and Montel, P. (eds),
"\OE uvres Compl\`etes" Gauthier-Villars, (1953)
\end{thebibliography}
\end{document}