#Program /home/stebla/programs/gap/Whiteheadv3.g to experiment
#with a GAP implementation of Whitehead's algebra.
#
#References
#[1] "A. N. Whitehead's Geometric Algebra", S. Blake, 19th February
#2005, http://stebla.pwp.blueyonder.co.uk/papers/Euclid.pdf
#[2] "Abstract Algebra", David S. Dummit and Richard M. Foote,
#Wiley International, 2004.
#History
#28th June 2005: Created. The idea here is to follow the way the 
#Whitehead.sci function package was created in Scilab. However, we 
#do not need to be as general as in the Scilab package because we 
#just need to make a wap function that works on the basis elements 
#represented as products of the reference points. Then, we can fill in 
#GAP's table of structure constants that define the algebra. In this
#GAP implementation the initial wap function that is used to fill in
#the SC table represents basis elements like a1a3a4 as [1,3,4].
#The function hyperplane form returns (say) +/-A1A3A4 as [sgn,[1,3,4]]
#where sgn is the sign. Similarly, the functions point_form and
#wap itself return such signed expressions in point form. This is 
#all we need to build the SC table.
#4th July 2005: Added elliptic, null and anti-linear polarities
#7th October 2006: Created a new version Whiteheadv2.g in order
#to experiment with defining the algebra over the field of rational
#functions or of Laurent polynomials in the indeterminate z=a1a2...an 
#instead of the complex numbers. The reason for doing this is because 
#the old (Whitehead.g) version assumed that a1a2..an=1 and this assumption 
#is not correct and gets us into all sorts of trouble whenever we use 
#functions f from which det(f) is not unity. It was necessary to delete
#all the polarities because they no longer work.
#12th October 2006: Rewrote the elliptic polarity to work with the new 
#understanding of Whitehead's algebra as a module over a polynomial ring.
#15th October 2006: Made a new version Whiteheadv3.g as a module over
#a polynomial ring, but using the basis a1,a2,...,a1a2...an. In Whiteheadv2.g
#the basis included 1 and excluded a1a2...an=z whilst in the new version
#the basis includes z=a1...an but excludes 1. The reason for this is that
#it may be the most natural version; work on the treatment of 
#reflections in Whitehead's algebra gives a single compact formula that 
#handles reflections for all point grades r=1 to r=n but breaks down for 
#the (ordinary) numbers of grade r=0. This suggests that the basis of the
#algebra should be over the grades r=1,...,n. 
wapv3:=function(F,R,n,x,y)
    #x and y are lists of the labels of reference points
    #F is the ground field and R is the polynomial ring and n is the 
    #number of reference points needed to span the space.
    local sgn,xy,z,x1,y2,x1x2,y1y2;
    z:=Indeterminate(F,1);
    if (Size(x)+Size(y) < n) then
        if Intersection(x,y)=[] then
            xy:=Concatenation(x,y);
            sgn:=SignPerm(MappingPermListList(AsSortedList(xy),xy));
            return [sgn*z^0,AsSortedList(xy)];
        else 
            return [Zero(R),[1..n]];
        fi;
    elif (Size(x)+Size(y) = n) then
        if Intersection(x,y)=[] then
            xy:=Concatenation(x,y);
            sgn:=SignPerm(PermList(xy));
            return [sgn*z^0,[1..n]];
        else 
            return [Zero(R),[1..n]];
        fi;
    elif (Size(x)+Size(y) > n) then
        if (Size(Intersection(x,y))=Size(x)+Size(y)-n) then
            xy:=Intersection(x,y);
            x1:=Difference(x,xy);
            y2:=Difference(y,xy);
            x1x2:=Concatenation(x1,xy);
            y1y2:=Concatenation(xy,y2);
            sgn:=SignPerm(MappingPermListList(x,x1x2));
            sgn:=sgn*SignPerm(MappingPermListList(y,y1y2));
            sgn:=sgn*SignPerm(PermList(Concatenation(x1x2,y2)));
            return [sgn*z,xy];
        else
            return [Zero(R),[1..n]];
        fi;
    fi;     
end;
#F:=CyclotomicField(4);
F:=Rationals;
R:=UnivariatePolynomialRing(F);
n:=4;#Grade of projective space
#The basis is a1,a2,...,a1a2...an. The pseudonumber a1a2...an
#should belong to the ring that acts on the vector space.
#We formalize this by setting z=a1a2..an as the indeterminate
#of a polynomial so that Whitehead's algebra is a module over
#a polynomial ring.
basis:=[];
for k in [1..n] do
    basis:=Concatenation(basis,Combinations([1..n],k));
od;
d:=Size(basis);
T:=EmptySCTable(d,Zero(R));#Create table for structure constants
for i in [1..d] do
    x:=basis[i];
    for j in [1..d] do
        y:=basis[j];
        xy:=wapv3(F,R,n,x,y);#xy[1]=sgn,xy[2]=element and sgn is a polynomial
        if xy[1]<>Zero(R) then #Fill in non-zero entries
            k:=Position(basis,xy[2]);
            SetEntrySCTable(T,i,j,[xy[1],k]);
        fi;
    od;
od;
#Make list ["a1","a2",...,"an"];
pts_names:=List([1..n],i->Concatenation("a",String(i)));
basis_names:=ShallowCopy(basis);
for i in [1..Length(basis)] do
    basis_names[i]:="";
    for j in basis[i] do
        basis_names[i]:=Concatenation(basis_names[i],pts_names[j]);
    od;
od;
WA:=AlgebraByStructureConstants(R,T,basis_names);
a:=Basis(WA);
#Here is an easy way to get the dual hyperplanes which uses the power
#of GAP to find the inverse, but it is computationally expensive.
#The dual element of a[1]..a[r] is A[1]..A[r] and the conventional way
#to get it is to write A[1]..A[r]=a[r+1]..a[n]/z where z=a[1]..a[n] because
#(a[1]..a[r])(A[1]..A[r])=1. In the first version of Whitehead.g, we used
#the function hyperplane_form to do this using GAP's permutation operations.
#However by setting x=a[1]..a[r] we can let GAP calculate the quotient 1/x
#which is defined so that (1/x)x=1. If we set the dual of x as y, then we
#want to solve xy=1 and y=(+/-)(1/x) and we can get the sign correctly
#by writing y=(1/x)(x(1/x)) because (x/(1/x)) is the sign and then,
#xy=(x(1/x))(x(1/x))=1.
A:=[];
for x in a do
    y:= (1/x)*(x*(1/x));
    Add(A,y);
od;
#
#Make the elliptic polarity
#It would have been nice to get the elliptic polarity I using 
#the following construct,
#polarity:=AlgebraHomomorphismByImages(WA,WA,a,A);
#in which I is an algebra/ring-module homomorphism. However, mathematicians 
#(see page 346-347 of [2]) define an R-module homomorphism (in our current
#notation) as I(rx+y)=rI(x)+I(y). In other words, they assume that the 
#ring elements r are scalars. However, in our case, the ring elements
#are polynomials r=p(z) where z=a[1]..a[n] and so I(rx+y)=I(r)I(x)+I(y).
#So, this function has to decompose the polynomial (assumed to be a 
#Laurent polynomial) and construct I(r).
I:=function(x)
    local p,Iz,z,Ix,c_i,Ip_i,i,j;
    p:=Coefficients(a,x);#Laurent polynomial coeffs
    Iz:=Product(A{[1..n]});z:=Indeterminate(F);
    Ix:=Zero(WA);
    for i in [1..Size(a)] do
        c_i:=CoefficientsOfLaurentPolynomial(p[i]);
        Ip_i:=Zero(WA);
        for j in [1..Size(c_i[1])] do
            #The coefficients c_i[1][j] are valued in the
            #ground field F, typically the Gaussian rationals.
            #GAP gives a no method found error on an attempt
            #to multiply (say) E(4)*a[1] whilst it does not 
            #complain about (say) 2*a[1]. So we have to
            #convert all our coefficients into Laurent
            #mononomials as c_i[1][j]*z^0 for the multiplication
            #to work.
            Ip_i:=Ip_i+(c_i[1][j]*z^0)*Iz^(c_i[2]+j-1);
        od;
        Ix:=Ix+Ip_i*A[i];
    od;
    return Ix;
end;