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

)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 ==> Join(AbelianGroup, Eltable(V(n,K),K)) 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.
 
    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))

)abbrev domain CLIF CliffordAlgebra
++ Author: Stephen M. Watt
++ Date Created: August 1988
++ Date Last Updated: May 17, 1991
++ Basic Operations: wholeRadix, fractRadix, wholeRagits, fractRagits
++ Related Domains: QuadraticForm, Quaternion, Complex
++ Also See:
++ AMS Classifications:
++ Keywords: clifford algebra, grassman algebra, spin algebra
++ Examples:
++ References:
++
++ Description:
++  CliffordAlgebra(n, K, Q) defines a vector space of dimension \spad{2**n}
++  over K, given a quadratic form Q on \spad{K**n}.
++
++  If \spad{e[i]}, \spad{1<=i<=n} is a basis for \spad{K**n} then
++     1, \spad{e[i]} (\spad{1<=i<=n}), \spad{e[i1]*e[i2]} 
++  (\spad{1<=i1<i2<=n}),...,\spad{e[1]*e[2]*..*e[n]}
++  is a basis for the Clifford Algebra.
++
++  The algebra is defined by the relations
++     \spad{e[i]*e[j] = -e[j]*e[i]}  (\spad{i \~~= j}),
++     \spad{e[i]*e[i] = Q(e[i])}
++
++  Examples of Clifford Algebras are: gaussians, quaternions, exterior 
++  algebras and spin algebras.
 
CliffordAlgebra(n, K, Q): T == Impl where
    n: PositiveInteger
    K: Field
    Q: QuadraticForm(n, K)
 
    PI ==> PositiveInteger
    NNI==> NonNegativeInteger
 
    T ==> Join(Ring, Algebra(K), VectorSpace(K)) with
        e: PI -> %
          ++ e(n) produces the appropriate unit element.
        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.
 
    Impl ==> add
        Qeelist :=  [Q unitVector(i::PositiveInteger) for i in 1..n]
        dim     :=  2**n
 
        Rep     := PrimitiveArray K
 
        New     ==> new(dim, 0$K)$Rep
 
        characteristic   == characteristic$K
        dimension()      == 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: Integer * x: % == (z := New; for i in 0..dim-1 repeat z.i := m*x.i; z)
        c: K * x: % == (z := New; for i in 0..dim-1 repeat z.i := c*x.i; z)
 
        0            == New
        1            == (z := New; z.0 := 1; z)
        coerce(m: Integer): % == (z := New; z.0 := m::K; z)
        coerce(c: K): % == (z := New; z.0 := c; z)
 
        e b ==
            b::NNI > n => error "No such basis element"
            iz := 2**((b-1)::NNI)
            z := New; z.iz := 1; z
 
        -- The ei*ej products could instead be precomputed in
        -- a (2**n)**2 multiplication table.
        addMonomProd(c1: K, b1: NNI, c2: K, b2: NNI, z: %): % ==
            c  := c1 * c2
            bz := b2
            for i in 0..n-1 | bit?(b1,i) repeat
                -- Apply rule  ei*ej = -ej*ei for i~=j
                k: NNI := 0
                for j in i+1..n-1 | bit?(b1, 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)::NNI
                else
                    bz:= bz + 2**i
            z.bz := z.bz + c
            z
 
        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 addMonomProd(x.ix,ix,y.iy,iy,z)
            z
 
        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: NNI := 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: NNI := 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]
 
        monomial(c, lb) ==
            r := canonMonom(c, lb)
            z := New
            z r.basel := r.coef
            z
        coefficient(z, lb) ==
            r := canonMonom(1, lb)
            r.coef = 0 => error "Cannot take coef of 0"
            z r.basel/r.coef
 
        Ex ==> OutputForm
 
        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
        coerce(x: %): Ex ==
            tl := [coerceMonom(x.i,i) for i in 0..dim-1 | not zero? x.i]
            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)

	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

