# File: QuatAlg.gap
# Description: GAP functions for computing fundamental
# regions for the action of the unit group of a
# quaternion algebra over the Rationals.
# Author: Assaf Wool
# email: assafwool@yahoo.com
# Date: July 2007
# Define the global quaternion algebra record
QuatGlob:=rec();
# Function to analyze D and initialize the quaternion object
findOrderType:=function(D)
local pr,i,len,p,q,r,ok;
if D<1 then
Print("ERROR: discriminant ",D," must be positive\n");
return 0;
fi;
pr:=FactorsInt(D);
# check for multiple primes
len:=Length(pr);
i:=1;
while i0 or len=0 then
Print("ERROR: discriminant has ",len," divisors\n");
return 0;
fi;
# find type, u, v
QuatGlob:=rec();
QuatGlob.disc:=D;
# find volume, 2 and 3 torsion, genus and number of generators
QuatGlob.vol:=1/6;
QuatGlob.n2:=1;
QuatGlob.n3:=1;
for p in pr do
QuatGlob.vol:=QuatGlob.vol*(p-1);
if p>2 then
QuatGlob.n2:=QuatGlob.n2*(1-Legendre(-1,p));
fi;
if p>3 then
QuatGlob.n3:=QuatGlob.n3*(1-Legendre(-3,p));
fi;
if p=2 then
QuatGlob.n3:=QuatGlob.n3*2;
fi;
od;
QuatGlob.twog:=2-QuatGlob.n2/2-QuatGlob.n3*2/3+
QuatGlob.vol;
QuatGlob.numgen:=QuatGlob.twog+QuatGlob.n2+QuatGlob.n3;
if QuatGlob.numgen>QuatGlob.twog then
QuatGlob.numgen:=QuatGlob.numgen-1;
fi;
ok:=0;
p:=pr[1];
q:=pr[2];
if p=2 and len=2 then
if Legendre(-1,q)=-1 then
QuatGlob.u:=1;
QuatGlob.v:=q;
QuatGlob.type:=1;
QuatGlob.den:=2;
ok:=1;
elif Legendre(-2,q)=-1 then
QuatGlob.u:=2;
QuatGlob.v:=q;
QuatGlob.type:=2;
QuatGlob.den:=2;
ok:=1;
fi;
fi;
if p<>2 and len=2 then
if Legendre(-1,p)=-1 and Legendre(-1,q)=-1 then
QuatGlob.u:=1;
QuatGlob.v:=D;
QuatGlob.type:=2;
QuatGlob.den:=2;
ok:=1;
elif Legendre(-p,q)=-1 and Legendre(q,p)=-1 then
QuatGlob.u:=p;
QuatGlob.v:=q;
QuatGlob.type:=2;
QuatGlob.den:=2;
ok:=1;
elif Legendre(-q,p)=-1 and Legendre(p,q)=-1 then
QuatGlob.u:=q;
QuatGlob.v:=p;
QuatGlob.type:=2;
QuatGlob.den:=2;
ok:=1;
fi;
fi;
if ok=1 then
QuatGlob.alg:=QuaternionAlgebra(Rationals,
-QuatGlob.u,QuatGlob.v);
QuatGlob.basis:=Basis(QuatGlob.alg);
return 1;
fi;
# general types
r:=5;
while true do
if Legendre(-1,r)=1 then
# this question takes care of both odd and even cases
if ForAll(pr,x->Legendre(x,r)=-1) then
QuatGlob.u:=D;
QuatGlob.v:=r;
QuatGlob.type:=3;
QuatGlob.den:=2*r;
QuatGlob.aux:=RootMod(-D,r);
QuatGlob.alg:=QuaternionAlgebra(Rationals,
-QuatGlob.u,QuatGlob.v);
QuatGlob.basis:=Basis(QuatGlob.alg);
return 1;
fi;
fi;
r:=NextPrimeInt(r);
od;
Print("ERROR: bad type, shouldn't be here\n");
return 0;
end;
# Function that finds all pairs (x,y) with norm of type x^2+Uy^2
# between min and max, appends them to the global lists and
# sorts by norm.
findComplexNorms:=function(U, max, min)
local x,y,n,order,tmpl;
# initialize u
# QuatGlob.u:=U;
# Initialize if starting from 0
if min=0 then
QuatGlob.x:=[];
QuatGlob.y:=[];
QuatGlob.norm:=[];
fi;
# Find the norms
x:=0;
while x*x<=max do
y:=0;
n:=x*x+U*y*y;
while n<=max do
if n>=min and (QuatGlob.type<>3 or RemInt(n,QuatGlob.v)=0) then
Append(QuatGlob.x,[x]);
Append(QuatGlob.y,[y]);
Append(QuatGlob.norm,[n]);
fi;
y:=y+1;
n:=x*x+U*y*y;
od;
x:=x+1;
od;
# Sort the lists according to norm
order:=Sortex(QuatGlob.norm);
tmpl:=ShallowCopy(QuatGlob.x);
QuatGlob.x:=Permuted(tmpl,order);
tmpl:=ShallowCopy(QuatGlob.y);
QuatGlob.y:=Permuted(tmpl,order);
end;
# Function that finds all quaternion elements in the
# algebra B(-U,V), in a maximal order defined by type.
# It uses the (x,y) pairs computed and sorted by norm.
findQuatElem:=function(V,type)
local pos,i,maxNm,tmpNm,a,b,c,d,tmpL,size,j,t;
# initialize v
# QuatGlob.v:=V;
# Initialize the element array
QuatGlob.elem:=[];
# make sure type is 1-3
if type<>1 and type<>2 and type<>3 then
Print("Error, type must be 1 or 2 or 3\n");
return;
fi;
maxNm:=QuatGlob.norm[Length(QuatGlob.norm)];
i:=1;
while QuatGlob.norm[i]*V0 and b<>0) then
j:=1;
while j<=size do
Add(tmpL, [tmpL[j][1],-tmpL[j][2],tmpL[j][3],tmpL[j][4]]);
j:=j+1;
od;
size:=size*2;
fi;
if (c<>0) then
j:=1;
while j<=size do
Add(tmpL, [tmpL[j][1],tmpL[j][2],-tmpL[j][3],tmpL[j][4]]);
j:=j+1;
od;
size:=size*2;
fi;
if (d<>0) then
j:=1;
while j<=size do
Add(tmpL, [tmpL[j][1],tmpL[j][2],tmpL[j][3],-tmpL[j][4]]);
j:=j+1;
od;
size:=size*2;
fi;
# filter elements according to order type
if (type=1 and RemInt(a,2)=RemInt(b,2) and
RemInt(a,2)=RemInt(c,2) and
RemInt(a,2)=RemInt(d,2)) or
(type=2 and RemInt(a,2)=RemInt(c,2) and
RemInt(b,2)=RemInt(d,2)) then
Append(QuatGlob.elem,tmpL/QuatGlob.den);
fi;
if type=3 then
for t in tmpL do
if RemInt(t[2],QuatGlob.v)=0 and
RemInt(t[2]-t[4],2)=0 and
RemInt(t[3]-QuatGlob.aux*(t[4]-t[2]),
QuatGlob.v)=0 and
RemInt(t[1]-t[3]-QuatGlob.aux*(t[2]-t[4]),
2*QuatGlob.v)=0 then
Add(QuatGlob.elem,t/QuatGlob.den);
fi;
od;
fi;
pos:=pos+1;
od;
i:=i+1;
od;
end;
# function that finds a list of non-redundant
# generators, using the element list. It finds
# all elements whose centers are not in any
# isometric circle of a previous element
findQuatGenerators:=function()
local i,len,n,Nm,Nm2,Nmg,first,Last,c,d,g,h,flag;
# initialize generators list
QuatGlob.gen:=[];
# compute all norms
len:=Length(QuatGlob.elem);
Nm:=[1..len];
Apply(Nm, i->QuatGlob.elem[i][3]^2+QuatGlob.u*QuatGlob.elem[i][4]^2);
# find first non zero norm, and find last index for each norm
first:=PositionNonZero(Nm);
Nm2:=0;
n:=0;
Last:=[];
for i in [first..len] do
if Nm[i]>Nm2 then
n:=i-1;
Nm2:=Nm[i];
fi;
Last[i]:=n;
od;
# go over elements and add to generators if
# no other element contains the center in its
# isometric circle.
n:=first;
while n<=len do
h:=QuatGlob.elem[n];
# first go over the current generators, they have the best
# chance of reducing the norm, this should improve
# speed
flag:=1;
for g in QuatGlob.gen do
Nmg:=g[3]*g[3]+QuatGlob.u*g[4]*g[4];
if Nmg=0 and d2[3]<0 then
return true;
fi;
if d1[3]<0 and d2[3]>=0 then
return false;
fi;
# case of y=0
if d1[3]=0 then
return (d1[2]>0 and (d2[3]<>0 or d2[2]<0));
fi;
if d2[3]=0 then
return (d2[2]<0 and (d1[3]<>0 or d1[2]>0));
fi;
# case of same y signs, no y=0
v1:=d1[2]/d1[3];
v2:=d2[2]/d2[3];
if v1>v2 then
return true;
fi;
return false;
end;
# Function that finds inverse pairs of generators
# and arranges them around the circle. Special
# care must be given for the case of U=1, where i
# is a special generator.
arrangeGenerators:=function()
local i,j,len,len2,g,dir,dir2,x,y,inv,map,gen,last,ng;
dir:=[];
i:=1;
len:=Length(QuatGlob.gen);
while i<=len do
g:=QuatGlob.gen[i];
x:=g[1]*g[3]-QuatGlob.u*g[2]*g[4];
y:=g[1]*g[4]+g[2]*g[3];
Add(dir,[i,x,y]);
i:=i+1;
od;
# sort along the circle
Sort(dir, cmpDir);
# Must do special treatment for u=1, remove
# generators with the same direction, replace
# generators by equivalent ones with inverse
# in the right half circle
if QuatGlob.u=1 then
dir2:=[];
if len>0 then
Add(dir2,dir[1]);
fi;
for i in [2..len] do
if cmpDir(dir[i-1],dir[i]) then
Add(dir2,dir[i]);
else
QuatGlob.gen[dir[i][1]]:=[1,1,0,0];
fi;
od;
dir:=dir2;
last:=-[1,dir[1][2],dir[1][3]];
i:=1;
g:=QuatGlob.gen[dir[i][1]];
while (not cmpDir(last,[1,dir[i][2],dir[i][3]])) do
x:=-g[1]*g[3]-g[2]*g[4];
y:=-g[1]*g[4]+g[2]*g[3];
# if inverse is in opposite half, multiply by i on the right
# and remember that u=1.
if ((not cmpDir([1,x,y],last)) or cmpDir([1,x,y],-last)) then
ng:=[-g[2],g[1],g[4],-g[3]];
if ng[1]<0 or (ng[1]=0 and ng[2]<0) then
QuatGlob.gen[dir[i][1]]:=-ng;
else
QuatGlob.gen[dir[i][1]]:=ng;
fi;
fi;
# special case of first/last element of order 2
g:=QuatGlob.gen[dir[i][1]];
if ((i=1 and g[1]=0) or
(not cmpDir([1,dir[i][2],dir[1][3]],last) and g[1]=0)) then
ng:=[-g[2],g[1],g[4],-g[3]];
if ng[1]<0 or (ng[1]=0 and ng[2]<0) then
QuatGlob.gen[dir[i][1]]:=-ng;
else
QuatGlob.gen[dir[i][1]]:=ng;
fi;
fi;
i:=i+1;
g:=QuatGlob.gen[dir[i][1]];
od;
while i<=Length(dir) do
QuatGlob.gen[dir[i][1]]:=[1,1,0,0];
i:=i+1;
od;
fi;
# find inverses
inv:=[];
for i in [1..len] do
inv[i]:=0;
od;
for i in [1..len] do
if inv[i]=0 then
if QuatGlob.gen[i][1]=0 then
inv[i]:=i;
else
for j in [i+1..len] do
if (QuatGlob.gen[i][1] = QuatGlob.gen[j][1] and
QuatGlob.gen[i][2] = -QuatGlob.gen[j][2] and
QuatGlob.gen[i][3] = -QuatGlob.gen[j][3] and
QuatGlob.gen[i][4] = -QuatGlob.gen[j][4]) then
inv[i]:=j;
inv[j]:=i;
break;
fi;
od;
fi;
fi;
od;
# create map for removing elements with no inverse
# and putting them in their order on the circle
len2:=Length(dir);
map:=[];
j:=1;
# make room for generator i if u=1
if QuatGlob.u=1 then
j:=2;
fi;
for i in [1..len2] do
if inv[dir[i][1]]<>0 then
map[dir[i][1]]:=j;
j:=j+1;
fi;
od;
# Put sorted generators and inverse information in the
# global object
QuatGlob.imap:=[];
gen:=[];
# put i in front if u=1
if QuatGlob.u=1 then
Add(QuatGlob.imap,1);
Add(gen,LinearCombination(QuatGlob.basis,[0,1,0,0]));
fi;
for i in [1..len2] do
if inv[dir[i][1]]<>0 then
Add(QuatGlob.imap,map[inv[dir[i][1]]]);
Add(gen,LinearCombination(QuatGlob.basis,QuatGlob.gen[dir[i][1]]));
fi;
od;
QuatGlob.gen:=gen;
end;
# Function for finding the relators
findRelations:=function()
local i,len,done,rel,g,curr,coef,ok,tor2,tor3,ng;
QuatGlob.rel:=[];
QuatGlob.grel:=[];
done:=[];
len:=Length(QuatGlob.gen);
# nothing to do if there are no generators
if len<1 then
return 0;
fi;
for i in [1..len] do
done[i]:=0;
od;
ng:=0;
for i in [1..len] do
if QuatGlob.imap[i]<>i then
ng:=ng+1/2;
else
ng:=ng+1;
fi;
od;
tor2:=0;
tor3:=0;
ok:=1;
# find relations by using the inverse and circle order
for i in [1..len] do
if done[i]=0 then
rel:=[i];
g:=QuatGlob.gen[i];
# special care of order 2 generators, add relation
coef:=Coefficients(QuatGlob.basis,g);
if coef[1]=0 then
Add(QuatGlob.rel,[i]);
Add(QuatGlob.grel,g);
tor2:=tor2+1;
fi;
done[i]:=1;
curr:=QuatGlob.imap[i]-1;
if curr<1 then
curr:=len;
fi;
while done[curr]=0 do
Add(rel,curr);
g:=g*QuatGlob.gen[curr];
done[curr]:=1;
curr:=QuatGlob.imap[curr]-1;
if curr<1 then
curr:=len;
fi;
od;
# add relation to global object, and check if it is a real
# relation.
Add(QuatGlob.rel,rel);
Add(QuatGlob.grel,g);
coef:=Coefficients(QuatGlob.basis,g);
if coef[1]=0 then
tor2:=tor2+1;
elif coef[1]=1/2 or coef[1]=-1/2 then
tor3:=tor3+1;
elif coef=[1,0,0,0] or coef=[-1,0,0,0] then
ng:=ng-1;
else
ok:=0;
fi;
fi;
od;
QuatGlob.tor2:=tor2;
QuatGlob.tor3:=tor3;
if tor2=0 and tor3=0 then
ng:=ng+1;
fi;
QuatGlob.ng:=ng;
if ok=1 and
[tor2,tor3,ng]<>[QuatGlob.n2,QuatGlob.n3,QuatGlob.numgen] then
ok:=0;
Print("WARN: Thought to be impossible (bug maybe?)\n",
[tor2,tor3,ng]," ",[QuatGlob.n2,QuatGlob.n3,QuatGlob.numgen],
"\n");
fi;
return ok;
end;
# function for finding a torsion element using w=i
# in case the group has no torsion.
findTorsionW:=function()
local w,i,done,len,g,flag,curr,coef,numl;
w:=QuatGlob.basis[2];
done:=[];
len:=Length(QuatGlob.gen);
for i in [1..len] do
done[i]:=0;
od;
i:=1;
flag:=true;
while flag and i<=len/2 do
if done[i]=0 then
g:=QuatGlob.gen[i];
done[i]:=1;
curr:=QuatGlob.imap[i];
if curr>len/2 then
g:=g*w;
curr:=curr-len/2;
fi;
curr:=curr-1;
if curr=0 then
g:=g*w;
curr:=curr+len/2;
fi;
while done[curr]=0 do
g:=g*QuatGlob.gen[curr];
done[curr]:=1;
curr:=QuatGlob.imap[curr];
if curr>len/2 then
g:=g*w;
curr:=curr-len/2;
fi;
curr:=curr-1;
if curr=0 then
g:=g*w;
curr:=curr+len/2;
fi;
od;
coef:=Coefficients(QuatGlob.basis,g);
if coef[1]=0 then
numl:=[NumeratorRat(coef[2]),NumeratorRat(coef[3]),NumeratorRat(coef[4])];
QuatGlob.wtor:=g/Gcd(numl);
return;
fi;
fi;
i:=i+1;
od;
Print("WARN: can't find w-torsion point\n");
end;
# function for reducing redundant generators
reduceGens:=function()
local i,len,inRel,j,coef,max,ind,norm,gc,
flag,inv,rlen,tmpl,rel,newRel,map,
tgen,tinv,trel,tgrel;
# go over relations, find where each generator is
# ignore order 2 generators relation
# Print(QuatGlob.rel,"\n",QuatGlob.imap,"\n");
len:=Length(QuatGlob.rel);
inRel:=[];
for i in [1..len] do
coef:=Coefficients(QuatGlob.basis,QuatGlob.grel[i]);
if not(Length(QuatGlob.rel[i])=1 and coef[1]=0) then
for j in QuatGlob.rel[i] do
inRel[j]:=i;
od;
fi;
od;
flag:=false;
while (not flag) do
# find relevant generator with
# maximal c^2+u*d^2 norm. Ignore generators
# whose inverse is in the same relation
max:=0;
ind:=0;
for i in [1..len] do
coef:=Coefficients(QuatGlob.basis,QuatGlob.grel[i]);
if (coef=[1,0,0,0] or coef=[-1,0,0,0]) then
for j in QuatGlob.rel[i] do
if (QuatGlob.imap[j]=j or
inRel[QuatGlob.imap[j]]<>i) then
gc:=Coefficients(QuatGlob.basis,QuatGlob.gen[j]);
norm:=gc[3]*gc[3]+QuatGlob.u*gc[4]*gc[4];
if norm>max then
max:=norm;
ind:=j;
fi;
fi;
od;
fi;
od;
if max>0 then
# update relations by substitution of gen[ind]
# Print("removing: ",ind, " ", QuatGlob.gen[ind]);
rel:=QuatGlob.rel[inRel[ind]];
i:=Position(rel,ind);
rlen:=Length(rel);
tmpl:=Concatenation(rel{[i+1..rlen]},rel{[1..i-1]});
# find the location of the inverse gen, and substitute
inv:=QuatGlob.imap[ind];
rel:=QuatGlob.rel[inRel[inv]];
i:=Position(rel,inv);
rlen:=Length(rel);
newRel:=Concatenation(rel{[1..i-1]},tmpl,rel{[i+1..rlen]});
# update inRel, and relations in QuatGlob
QuatGlob.rel[inRel[ind]]:=[];
QuatGlob.rel[inRel[inv]]:=newRel;
# Print(" ", newRel);
inRel[ind]:=0;
for i in tmpl do
inRel[i]:=inRel[inv];
od;
inRel[inv]:=0;
else
flag:=true;
fi;
od;
# find removed generators
j:=1;
map:=[];
len:=Length(QuatGlob.gen);
for i in [1..len] do
if inRel[i]<>0 then
Add(map,j);
j:=j+1;
else
Add(map,0);
fi;
od;
# update generators and inverse
tgen:=[];
tinv:=[];
for i in [1..len] do
if inRel[i]<>0 then
Add(tgen,QuatGlob.gen[i]);
Add(tinv,map[QuatGlob.imap[i]]);
fi;
od;
QuatGlob.gen:=tgen;
QuatGlob.imap:=tinv;
# update relations
trel:=[];
tgrel:=[];
i:=1;
len:=Length(QuatGlob.rel);
while i<=len do
if Length(QuatGlob.rel[i])>0 then
Add(tgrel,QuatGlob.grel[i]);
Apply(QuatGlob.rel[i],j->map[j]);
Add(trel,QuatGlob.rel[i]);
fi;
i:=i+1;
od;
QuatGlob.grel:=tgrel;
QuatGlob.rel:=trel;
end;
# find and print all vertices
findVertices:=function()
local coef,len,i,j,g,cg,w,x,y,first,co,
dir,tc,ind,pos,rel,rlen,norm,allDir,
ok,num,den;
len:=Length(QuatGlob.rel);
norm:=1;
allDir:=[];
for i in [1..len] do
g:=QuatGlob.grel[i];
rel:=QuatGlob.rel[i];
coef:=Coefficients(QuatGlob.basis,g);
if coef=[1,0,0,0] or coef=[-1,0,0,0] then
if len > 1 then
Print("ERROR: more than one surface relation\n");
return;
fi;
if not IsBound(QuatGlob.wtor) then
Print("ERROR: no w-torsion point found\n");
return;
fi;
# init g and find the starting point for the relation orbit
g:=QuatGlob.wtor;
coef:=Coefficients(QuatGlob.basis,g);
norm:=coef[1]^2+QuatGlob.u*coef[2]^2-
QuatGlob.v*(coef[3]^2+QuatGlob.u*coef[4]^2);
# be careful about the direction of the w-torsion
# point, sign of the real coefficient has to be considered
if (coef[2]>0 and coef[2]^2*QuatGlob.u>norm) or
(coef[2]<0 and coef[2]^2*QuatGlob.u0 and cg[2]^2*QuatGlob.u>norm) or
(cg[2]<0 and cg[2]^2*QuatGlob.u0 then
Print("-S(",norm,"))");
else
Print("+S(",norm,"))");
fi;
Print("/", "S(",QuatGlob.u*QuatGlob.v,") x ");
Print("(",-QuatGlob.u*cg[4]/num*den);
if cg[3]>=0 then
Print("+");
fi;
Print(cg[3]/num*den,"S(-",QuatGlob.u,"))\n");
# print the relation coefficients too
Print(j,": rel " , cg[1], " ", cg[2], " ", cg[3], " ", cg[4],"\n");
fi;
# shift relation cyclically
g:=QuatGlob.gen[QuatGlob.imap[j]]*g*QuatGlob.gen[j];
od;
elif coef[1]=1/2 or coef[1]=-1/2 then
for j in rel do
cg:=Coefficients(QuatGlob.basis,g);
if (cg[2]>0 and 4*cg[2]^2*QuatGlob.u>3) or
(cg[2]<0 and 4*cg[2]^2*QuatGlob.u<3) then
allDir[j]:=[1,-QuatGlob.u*cg[4],cg[3]];
else
allDir[j]:=[1,QuatGlob.u*cg[4],-cg[3]];
fi;
Print(j,": (",2*cg[2],"S(",QuatGlob.u,")");
if cg[2]>0 then
Print("-S(3))");
else
Print("+S(3))");
fi;
Print("/(",2*(cg[3]*cg[3]+QuatGlob.u*cg[4]*cg[4]),
"S(",QuatGlob.u*QuatGlob.v,") x ");
Print("(",-QuatGlob.u*cg[4]);
if cg[3]>=0 then
Print("+");
fi;
Print(cg[3],"S(-",QuatGlob.u,"))\n");
# print the relation coefficients too
Print(j,": rel " , cg[1], " ", cg[2], " ", cg[3], " ", cg[4],"\n");
# shift relation cyclically
g:=QuatGlob.gen[QuatGlob.imap[j]]*g*QuatGlob.gen[j];
od;
else
Print("ERROR: non-torsion relation ",rel,"\n");
return;
fi;
od;
# verify that vertices are ordered as they should
ok:=true;
len:=Length(QuatGlob.gen);
for i in [2..len-2] do
if (not cmpDir(allDir[i],allDir[i+1])) then
ok:=false;
break;
fi;
od;
if not ( (cmpDir(allDir[1],allDir[2]) and
cmpDir(allDir[len],allDir[1]) ) or
(cmpDir(allDir[1],allDir[2]) and
cmpDir(allDir[len-1],allDir[len]) ) or
(cmpDir(allDir[len],allDir[1]) and
cmpDir(allDir[len-1],allDir[len]) ) ) then
ok:=false;
fi;
if (not ok) then
Print("ERROR: vertices are not sorted along the circle\n");
fi;
end;
# Runner function, take D and get a full object with
# generators and fundamental region.
createQuatObj:=function(D)
local ok,max,min;
ok:=findOrderType(D);
if ok=0 then
return 0;
fi;
max:=QuatGlob.disc*QuatGlob.disc;
min:=0;
ok:=0;
while ok=0 do
findComplexNorms(QuatGlob.u, max, min);
findQuatElem(QuatGlob.v,QuatGlob.type);
findQuatGenerators();
arrangeGenerators();
ok:=findRelations();
if ok=0 then
min:=max+1;
max:=2*max;
Print("not enough generators, increase max norm to ",max,"\n");
fi;
od;
# use w=i to find torsion element in group with w
# if group has no torsion.
if QuatGlob.n2=0 and QuatGlob.n3=0 then
findTorsionW();
fi;
# find redundant generators until all relations are eliminated
# except torsion or surface relations
reduceGens();
# find and print vertices of region
findVertices();
return 1;
end;
Text file Source (historic): geocities.com/assafwool/Quat
geocities.com/assafwool(to report bad content: archivehelp @ gmail)
|
|
|
|
|