build open-axiom
--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.
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 ==
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 := members(ans)
ansC:% := 0
for pj in Pn repeat
cj:= first ansP
ansP := rest ansP
ansC := ansC + cj*monomial(1,pj)
ansC