Quelle  bianchi.gi   Sprache: unbekannt

#Throughout we try to denote Sqrt(-d) by S

############################################################
############################################################
InstallGlobalFunction(HAP_BianchiFundamentalRectangle, #####
function(OQ)                                           #####
local s1, s2, t1, t2;                       
#The fundamental rectangle is [s1,t1*S]x[s2,t2*S].
#The numbers s1,t1,s2,t2 are returned. We also use the 
#subregion Q=[0,t1*S]x[0,t2*S].

if not -OQ!.bianchiInteger mod 4 = 3  then
s1:=-1/2;t1:=1/2;s2:=-1/2;t2:=1/2;
else
s1:=-1/2;t1:=1/2;s2:=-1/4;t2:=1/4;
fi;
return [s1,t1,s2,t2];                                  #####
end);                                                  #####
############################################################ 
############################################################


############################################################
############################################################
InstallGlobalFunction(UnimodularPairStandardForm,      #####
function(x)                                            #####
local y;                                    
#The unimodular pairs [a,b] and [-a,-b] are equivalent.
#One of these pairs is returned.
y:=1*x;
y[1]:=-y[1];
y[2]:=-y[2];
return Maximum(x,y);;                                  #####
end);                                                  #####
############################################################
############################################################


####################################################
####################################################
InstallGlobalFunction(QuadraticIntegersByNorm, #####
function(OQ,N)      
local BI, NW, L, x,y,z;                      
#Returns all quadratic integers a with 
#-N <= Norm(OQ,a) <= N where OQ is a ring of integers

if N<0 then return []; fi;

BI:=OQ!.bianchiInteger;

if IsBound(OQ!.normomega) then NW:=OQ!.normomega; 
else
    if -BI mod 4 =3 then NW:=(1+AbsInt(BI))/4; OQ!.normomega:=NW;
    else NW:=AbsInt(BI); OQ!.normomega:=NW;
    fi;
fi;

L:=[]; #To be a list of integers of appropriate norm.

if -BI mod 4 =3 then
for x in [-(1+Int(Sqrt(1.0*N)))..1+Int(Sqrt(1.0*N))] do
for y in [-(1+Int(Sqrt(1.0*N/NW)))..1+Int(Sqrt(1.0*N/NW))] do
z:=QuadraticNumber((2*x+y)/2,y/2,BI);  
if HAPNorm(OQ,z)<=N then Add(L,z); fi;
od;
od;
else
for x in [-(1+Int(Sqrt(1.0*N)))..1+Int(Sqrt(1.0*N))] do
for y in [-(1+Int(Sqrt(1.0*N/NW)))..1+Int(Sqrt(1.0*N/NW))] do
z:=QuadraticNumber(x,y,BI);
if HAPNorm(OQ,z)<=N then Add(L,z); fi;
od;
od;
fi;

SortBy(L,x->HAPNorm(OQ,x));

return L;                                    #####
end);                                        #####
##################################################
##################################################


########################################################
########################################################
InstallGlobalFunction(HAP_UnimodularComplements,   #####
function(OQ,a)                                     #####
local Q, nrma,BI, rat, irrat, N, L, U, b, x,s1, t1, s2, t2;
#Returns the list of unimodular pairs (a,b) in OQxOQ that 
#intersect with the fundamental rectangle Q. The pair (a,b) 
#is returned as [a,b,Norm(a)]. Here we think of (a,b) as 
#being represented by a circular disk with centre b/a and radius |a| 
#in the complex plane.


BI:=OQ!.bianchiInteger;

Q:=HAPQuadratic(a);
nrma:=HAPNorm(OQ,a);

L:=QuadraticIntegersByNorm(OQ,(1-BI)*nrma);   

L:=Filtered(L,x->not x=0);
U:=[];
for b in L do
if IsQQUnimodularPair(OQ,[a,b,nrma]) then

if Norm(QuadraticIdeal(OQ,[QuadraticToCyclotomic(a),QuadraticToCyclotomic(b)])) = 1 
then Add(U,UnimodularPairStandardForm([a,b,nrma])); fi;

fi;

#fi;
od;
#Add(U,[1,0,1]);      #Added July 2025
return SSortedList(U);                       #####
end);                                        #####
##################################################
##################################################


#######################################################
#######################################################
InstallGlobalFunction(HAP_SqrtInequality,         #####
function(a,b,c);                                  #####
#Inputs three positive rational numbers and       
#returns Sqrt(a) + Sqrt(b) <= Sqrt(c).
if a+b>c then return false; fi;
return  4*a*b <= (c-a-b)^2;                       #####
end);                                             #####
#######################################################
#######################################################


#######################################################
#######################################################
InstallGlobalFunction(HAP_SqrtStrictInequality,   #####
function(a,b,c);                                  #####
#Inputs three positive rational numbers and      
#returns Sqrt(a) + Sqrt(b) <= Sqrt(c).
if a+b>=c then return false; fi;
return  4*a*b < (c-a-b)^2;                        #####
end);                                             #####
#######################################################
#######################################################


##########################################################
##########################################################
InstallGlobalFunction(HAP_IsRedundantUnimodularPair, #####
function(OQ,L,x)                                     #####
local y,l,m,lm,nrm,a,b,ba,rm,ra,nrmm,nrma;        
#Inputs a ring of integers OQ, a list L of unimodular 
#pairs in OQ, and a unimodular pair x. Returns true if 
#x is 'redundant' and false otherwise. We think of a
#unimodular pair [a,b] as a sphere of radius 1/|a| and centre
#b/a in the complex plane. So we assume not a=0. A sphere is
#redundant in this function if it lies underneath some other 
#sphere.

l:=x[2];
m:=x[1];
nrmm:=x[3];

for y in L do
a:=y[1];
b:=y[2];
nrma:=y[3];

nrm:=HAPNorm(OQ,l*a - b*m)/(nrmm*nrma);
if not nrm > 1/nrma then #so circle y might contain circle x
if HAP_SqrtInequality(nrm,1/nrmm,1/nrma) then return true; fi;
fi;
od;

return false;                                #####
end);                                        #####
##################################################
##################################################


######################################################
######################################################
InstallGlobalFunction(UnimodularPairCoordinates, #####
function(OQ,x)                                   #####
local a,b,Q,ax,ay,bx,by, D,X,Y,BI;          
#inputs a unimodular pair [a,b,1/r^2] and returns [x1,x2,r^2] 
#where r is the radius and [x1,x2*S] the centre of the 
#corresponding hemisphere. Here S:=Sqrt(BianchiInteger)

if Length(x)>3 then return [x[4],x[5],x[6]]; fi;

a:=x[2];
b:=x[1];
BI:=OQ!.bianchiInteger;

Q:=HAPQuadratic(2*a);
ax:=Q.a/(2*Q.d);
ay:=Q.b/(2*Q.d);
Q:=HAPQuadratic(2*b);
bx:=Q.a/(2*Q.d);
by:=Q.b/(2*Q.d);

D:=bx^2 - by^2*BI;
X:=ax*bx-ay*by*BI;
Y:=ay*bx-ax*by;

D:=D^-1;
x[4]:=X*D;
x[5]:=Y*D;
x[6]:=1/x[3];

return  [x[4],x[5],x[6]];
                                             #####
end);                                        #####
##################################################
##################################################


##################################################
##################################################
InstallGlobalFunction(UnimodularPairs,       #####
function(arg)                                #####
local OQ,N,bool,BOOL,L,LL,a,b,fn,U,UU,u,w,   
K,i,j,nrm,NRM;
#Returns a list of pairs [a,b] with Norm(a)<=N and (a,b)
#a unimodular pair intersecting the closed rectangle Q=
#[0,t1*S]x[0,t2*S]. With arg[3]=false no redundancies are 
#eliminated. With arg[4]=LL a list of unimodular pairs up 
#to a given norm -- some work is avoided.

OQ:=arg[1];
N:=arg[2];
if Length(arg)>2 then bool:=arg[3]; else bool:=true; fi;
#bool=true means that some obvious redundancies will be removed.
if Length(arg)>3 then LL:=List(arg[4],x->1*x); else LL:=[]; fi;
#LL is a list of all unimodular pairs of a neighbourhood of the rectangle 
#up to a certain norm.
if Length(arg)>4 then BOOL:=arg[5]; else BOOL:=false; fi;
#BOOL=true means a slightly larger neighbourhood of the fundamental recangle is considered.


L:=QuadraticIntegersByNorm(OQ,N);
L:=Filtered(L,x->not x=0);
if Length(LL)>0 then
NRM:=Maximum(List(LL,x->x[3]));
L:=Filtered(L,x->HAPNorm(OQ,x)>NRM);
else NRM:=0;
fi;

for a in L do
U:=HAP_UnimodularComplements(OQ,a);

if bool then
UU:=[];
for u in U do
  if not HAP_IsRedundantUnimodularPair(OQ,LL,u) then
  Add(UU,u);
  fi;
od;
Append(LL,UU);
else
Append(LL,U);
fi;
od;

#LL:=SSortedList(LL);
LL:=Classify(LL,x->x);;
LL:=List(LL,x->x[1]);
UU:= Filtered(LL,x->x[3]>NRM);
Append(LL,UU);

if bool then
for i in Reversed([2..Length(LL)]) do
nrm:=LL[i][3];
for j in Reversed([1..i-1]) do
if LL[j][3] < nrm then
if HAP_IsRedundantUnimodularPair(OQ,LL{[1..j]},LL[i]) then
LL[i]:=0;
fi;
break;
fi;
od;
od;
fi;

LL:=Filtered(LL,x->not x=0);
LL:=List(LL,x->[x[1],x[2],HAPNorm(OQ,x[1])]);
#LL:=SSortedList(LL);
LL:=Classify(LL,x->x);;
LL:=List(LL,x->x[1]);

SortBy(LL,x->x[3]);
#List(LL,x->UnimodularPairCoordinates(OQ,x));
return LL;                                   #####
end);                                        #####
##################################################
##################################################


#######################################################
#######################################################
InstallGlobalFunction(UnimodularIntersectingLine, #####
function(OQ,a,b)                                  #####
local ca,cb,u,v,w,S,ra2,rb2,u2,lambda;
#Inputs two unimodular pairs a, b representing circles
#in the complex plane. It is assumed that the circles 
#intersect. A pair [v,w] is returned where w is the 
#direction vector of the radical axis and v is the point 
#where the radical acis intersects the line joining 
#the centres of a and b. 

S:=QuadraticNumber(0,1,-OQ!.bianchiInteger);

ca:=UnimodularPairCoordinates(OQ,a);
ca[2]:=ca[2]*S;
cb:=UnimodularPairCoordinates(OQ,b);
cb[2]:=cb[2]*S;
u:=cb-ca;
u:=u{[1,2]};
w:=[u[2],-u[1]];

if u=[0,0] then return [ca{[1,2]},w]; fi;

ra2:=ca[3];
rb2:=cb[3];
#u2:=Norm(OQ,u[1])+Norm(OQ,u[2]);
u2:=u[1]^2+u[2]^2;

ca:=ca{[1,2]};
lambda:=(ra2-rb2+u2)/(2*u2);

v:=ca+lambda*u;
return [v,w];                                     #####
end);                                             #####
#######################################################
#######################################################


#########################################################
#########################################################
InstallGlobalFunction(UnimodularIntersectingPoint,  #####
function(OQ,a,b,c,Lab,Lac)                          #####
local v,w,vv,ww,D,lambda;
#Inputs three unimodular pairs a,b,c thought of as
#circles in the complex plane. It is assumed that the
#circles intersect pairwise. The radical centre v
#is returned.

#Lab:=UnimodularIntersectingLine(OQ,a,b);
#Lac:=UnimodularIntersectingLine(OQ,a,c);

v:=Lab[1];
w:=Lab[2];
vv:=Lac[1];
ww:=Lac[2];
D:=w[1]*(-ww[2])+w[2]*ww[1];
if D=0 then
return fail; fi;
D:=D^-1;
lambda:=D*( -ww[2]*(vv[1]-v[1]) + ww[1]*(vv[2]-v[2]) );

return v+lambda*w;                                #####
end);                                             #####
#######################################################
#######################################################


############################################################
############################################################
InstallGlobalFunction(HAP_HeightOfPointOnSphere,       #####
function(OQ,z,p)                                       #####
local c,S;                                          
#Inputs the real and imaginary coordinates [x,y] of a point
#z=x+i*y in the complex plain and unimodular pair p. 
#The pair p determines a hemisphere centred at some point 
#on the complex plane. The square h^2 of the height h of 
#the lift of z to the hemisphere is returned. If z does 
#not lie under the hemisphere then no lift exists and a 
#negative square height will be returned. 

S:=QuadraticNumber(0,1,-OQ!.bianchiInteger);

c:=UnimodularPairCoordinates(OQ,p);
S:=c[3]-(z[1]-c[1])^2-(z[2]-S*c[2])^2;
return S;                                             #####
end);                                                 #####
###########################################################
###########################################################


################################################################
################################################################
InstallGlobalFunction(HAP_AreIntersectingUnimodularPairs,  #####
function(OQ,a,b)                                           #####
local A, B, BI;                                        
#Inputs a ring of integers OQ of an imaginary quadratic field
#together with two unimodular pairs.
#It returns true if the two circles in the commplex plane
#corresponding to the unimodular pairs intersect.  

BI:=AbsInt(OQ!.bianchiInteger);

A:=UnimodularPairCoordinates(OQ,a);
B:=UnimodularPairCoordinates(OQ,b);
return not HAP_SqrtStrictInequality(A[3],B[3], (A[1]-B[1])^2 + BI*(A[2]-B[2])^2 );
                                                       #####
end);                                                  #####
############################################################
############################################################


################################################################
################################################################
InstallGlobalFunction(HAP_AreStrictlyIntersectingUnimodularPairs,  #####
function(OQ,a,b)                                           #####
local A, B, BI;
#Inputs a ring of integers OQ of an imaginary quadratic field
#together with two unimodular pairs.
#It returns true if the two circles in the commplex plane
#corresponding to the unimodular pairs intersect.

BI:=AbsInt(OQ!.bianchiInteger);

A:=UnimodularPairCoordinates(OQ,a);
B:=UnimodularPairCoordinates(OQ,b);
return not HAP_SqrtInequality(A[3],B[3], (A[1]-B[1])^2 + BI*(A[2]-B[2])^2 );
                                                       #####
end);                                                  #####
############################################################
############################################################



#################################################################
#################################################################
InstallGlobalFunction(HAP_Are3IntersectingUnimodularPairs,  #####
function(arg)                                               #####
local OQ, a,b,c, Lab, Lac, Lbc, A, B, C, v, w;                
#Inputs a ring of integers OQ of an imaginary quadratic field
#together with three unimodular pairs that intersect pairwise.
#It returns true if the three circles in the commplex plane
#corresponding to the unimodular pairs have common region 
#of overlap. The common intersections may be a point. 
#Otherwise false is returned.

OQ:=arg[1];
a:=arg[2];
b:=arg[3];
c:=arg[4];
Lab:=arg[5];
Lac:=arg[6];
Lbc:=arg[7];

###########################
A:=UnimodularPairCoordinates(OQ,a);
B:=UnimodularPairCoordinates(OQ,b);
C:=UnimodularPairCoordinates(OQ,c);

if HAP_HeightOfPointOnSphere(OQ,Lab[1],c)>=0 then return [true]; fi;

if HAP_HeightOfPointOnSphere(OQ,Lac[1],b)>=0 then return [true]; fi;

v:=Lab[2]{[2,1]};
v[1]:=-v[1];
if v*Lac[2]=0 then return [true]; fi;    #circles are colinear

#v:=UnimodularIntersectingLine(OQ,c,b);
if HAP_HeightOfPointOnSphere(OQ,Lbc[1],a)>=0 then return [true]; fi;

v:=UnimodularIntersectingPoint(OQ,a,b,c,Lab,Lac);

if HAP_HeightOfPointOnSphere(OQ,v,a)<0 then return [false]; fi;
if HAP_HeightOfPointOnSphere(OQ,v,b)<0 then return [false]; fi;
if HAP_HeightOfPointOnSphere(OQ,v,c)<0 then return [false]; fi;

return [true];                               #####
end);                                        #####
##################################################
##################################################


##################################################
##################################################
InstallGlobalFunction(HAP_PrintFloat,        #####
function(xx)                                 #####
local x,s,m, p, l, sgn;                      
x:=xx;
sgn:=SignFloat(x);
x:=AbsoluteValue(x);

s:=SplitString(String(x),'e');;
if Length(s)=1 then
s:=SplitString(String(xx),'e');;
return s[1]{[1..Minimum(5,Length(s[1]))]}; fi;

p:=EvalString(s[2]);;
l:=Length(s[1]);
m:=s[1]{[2..l]};
if p<0 then
s:=Concatenation(ListWithIdenticalEntries(-p,"0"));
s:=Concatenation(".",s,m{[1..l-1]});
s:=s{[1..5]};
if sgn<0 then s:=Concatenation("-",s); fi;
return s;
fi;
m:=EvalString(m);
s:=10^p*m;
s:=Concatenation(String(s),"00000");
s:=Concatenation(s{[1..p]},".",s{[p+1..Minimum(l,p+5)]});
if sgn<0 then s:=Concatenation("-",s); fi;
return s;                                    #####
end);                                        #####
##################################################
##################################################


###################################################
###################################################
InstallGlobalFunction(DisplayUnimodularPairs, #####
function(arg)                                #####
local OQ,L, C, file, i, x, w, BI,ABI, S1,S,SS,s1,t1,s2,t2;

OQ:=arg[1];
L:=arg[2];
if Length(arg)=3 then C:=arg[3]; 
else
C:=[];
fi;


BI:=OQ!.bianchiInteger;
ABI:=AbsInt(BI);
S:=Sqrt(1.0*ABI);
S1:=Sqrt(1.0*ABI);
if -BI mod 4 =3 then S:=0.25*S; else S:=0.5*S; fi;
SS:=HAP_PrintFloat(S);
SS:=SS{[1..Minimum(Length(SS),5)]};

C:=List(C,x->[x[1], x[2]*(1/QuadraticNumber(0,1,ABI))]);

C:=List(C,x->[x[1]!.rational,x[2]!.rational]);

C:=List(C,x->[1.0*x[1],S1*x[2]]);

file:="tmpp.asy";
PrintTo(file,"size(0,200);\n");
AppendTo(file,"defaultpen(0.2);\n");
AppendTo(file,"pen colour=blue+opacity(0.50);\n");
AppendTo(file,"real S=");
AppendTo(file,SS);
AppendTo(file,";\n");

s1:=HAP_PrintFloat(-0.5);t1:=HAP_PrintFloat(0.5);s2:=HAP_PrintFloat(-S);t2:=HAP_PrintFloat(S);

AppendTo(file,"path g=(-0.5,-S)--(0.5,-S)--(0.5,S)--(-0.5,S)--cycle;\n");

AppendTo(file,"picture fd;\n");

S:=Sqrt(1.0*ABI);
for x in L do
w:=UnimodularPairCoordinates(OQ,x);
w:=1.0*w;
AppendTo(file,"real r=");
AppendTo(file,"sqrt(");
AppendTo(file,HAP_PrintFloat(w[3]));
AppendTo(file,");\n");
AppendTo(file,"pair z=(");
AppendTo(file,HAP_PrintFloat(w[1]));
AppendTo(file,",");
AppendTo(file,HAP_PrintFloat(w[2]*S));
AppendTo(file,");\n");
AppendTo(file,"path c=circle(z,r);\n");
AppendTo(file,"fill(fd,c,colour);");
AppendTo(file,"draw(fd,c);\n");
od;

#AppendTo(file,"clip(fd,g);\n");
AppendTo(file,"draw(fd,g,white+linewidth(.3));\n");

for x in C do
AppendTo(file,"real r=0.01");
AppendTo(file,";\n");
AppendTo(file,"pair z=(");
AppendTo(file,HAP_PrintFloat(x[1]));
AppendTo(file,",");
AppendTo(file,HAP_PrintFloat(x[2]));
AppendTo(file,");\n");
AppendTo(file,"path c=circle(z,r);\n");
AppendTo(file,"fill(fd,c,red);");
AppendTo(file,"draw(fd,c);\n");
od;

AppendTo(file,"add(fd);\n");


Exec("asy -f pdf tmpp.asy; evince tmpp.pdf");     #####
end);                                              #####
#######################################################
#######################################################



#####################################################
#####################################################
InstallGlobalFunction(Display3DUnimodularPairs, #####
function(OQ,LL)                                 #####
local  file, x, w, BI,ABI, S,SS, L,R,T,B,sz,V;

V:=HAP_VertexHeights(OQ,LL)[2];
V:=V!.Heights;

V:=Filtered(V,x->x>0);

V:=List(V,x->x!.rational);

V:=1.0*Minimum(V);
V:=Sqrt(V);
V:=Minimum(0.05,V/2);


BI:=OQ!.bianchiInteger;
ABI:=AbsInt(BI);
S:=Sqrt(1.0*ABI);
if -BI mod 4 = 3 then S:=0.25*S; else S:=0.5*S; fi;
SS:=HAP_PrintFloat(S);
SS:=SS{[1..Minimum(Length(SS),5)]};
sz:=SizeScreen();
SizeScreen([200,24]);
file:="tmpp.asy";
PrintTo(file,"import solids;\n import graph3;\n import three;\n size(1000);\n");
AppendTo(file,"currentprojection=orthographic(5,4,10);\n");
AppendTo(file,"currentlight=Headlamp;\n");
AppendTo(file,"nslice=4*nslice;\n");
AppendTo(file,"revolution b=sphere(O,1);\n");

L:=[];R:=[];T:=[];B:=[];
S:=Sqrt(1.0*ABI);
for x in LL do
w:=UnimodularPairCoordinates(OQ,x);
w:=1.0*w;
w[2]:=w[2]*S;
Add(L,w[1]-Sqrt(w[3]));
Add(R,w[1]+Sqrt(w[3]));
Add(T,w[2]+Sqrt(w[3]));
Add(B,w[2]-Sqrt(w[3]));
AppendTo(file,"real r=");
AppendTo(file,"sqrt(");
AppendTo(file,HAP_PrintFloat(w[3]));
AppendTo(file,");\n");
AppendTo(file,"draw(shift(",HAP_PrintFloat(w[1]),",",HAP_PrintFloat(w[2]),",0)*scale3(r)*surface(b),blue+opacity(1));\n");
od;


AppendTo(file,"real S=");
AppendTo(file,SS);
AppendTo(file,";\n");
AppendTo(file,"real V=");
AppendTo(file,V);
AppendTo(file,";\n");

#if -BI mod 4 = 3 then
AppendTo(file,"path3 g=(-0.5,-S,1)--(0.5,-S,1)--(0.5,S,1)--(-0.5,S,1)--cycle;\n");
#else
#AppendTo(file,"path3 g=(0,0,1)--(1,0,1)--(1,S,1)--(0,S,1)--cycle;\n");
#AppendTo(file,"path3 gg=(0,0,V)--(1,0,V)--(1,S,V)--(0,S,V)--cycle;\n");
#fi;

AppendTo(file,"draw(g,white+linewidth(0.3));");

AppendTo(file,"real L=");
AppendTo(file,Minimum(L));
AppendTo(file,";\n");

AppendTo(file,"real R=");
AppendTo(file,Maximum(R));
AppendTo(file,";\n");

AppendTo(file,"real T=");
AppendTo(file,Maximum(T));
AppendTo(file,";\n");

AppendTo(file,"real B=");
AppendTo(file,Minimum(B));
AppendTo(file,";\n");


AppendTo(file,"triple p0 =(L,T,V), p1=(R,T,V), p2 = (R,B,V), p3 = (L,B,V);\n");
AppendTo(file,"triple q0 =(L,T,-1), q1=(R,T,-1), q2 = (R,B,-1), q3 = (L,B,-1);\n");

AppendTo(file,"skeleton s;\n");
AppendTo(file,"b.transverse(s,reltime(b.g,0.5),P=currentprojection);\n");
AppendTo(file,"draw(surface(p0--p1--p2--p3--cycle),white);\n");
AppendTo(file,"draw(surface(q0--q1--q2--q3--cycle),white);\n");
AppendTo(file,"draw(surface(p0--p1--q1--q0--cycle),white);\n");
AppendTo(file,"draw(surface(q1--p1--p2--q2--cycle),white);\n");
AppendTo(file,"draw(surface(q3--q2--p2--p3--cycle),white);\n");
AppendTo(file,"draw(surface(p0--q0--q3--p3--cycle),white);\n");

SizeScreen(sz);

Exec("asy -V tmpp.asy");                           #####
end);                                              #####
########################################################
########################################################


###################################################
###################################################
InstallGlobalFunction(IsQUnimodularPair,     #####
function(arg)                                 #####
local OQ,x,BI,ABI,L,s1,t1,s2,t2,T,S;                
#Inputs a unimodular pair x=[a,b] and returns true
#if the associated hemisphere touches the 
#fundamental rectangle [s1,t1*S]x[s2,t2*S] in the 
#complex plane.


OQ:=arg[1];
BI:=OQ!.bianchiInteger;
ABI:=AbsInt(BI);
x:=arg[2];
L:=HAP_BianchiFundamentalRectangle(OQ);
s1:=L[1];t1:=L[2];; s2:=L[3]; t2:=L[4];

T:=UnimodularPairCoordinates(OQ,x);

if T[1]>=0 and
(not  HAP_SqrtStrictInequality(t1^2,T[3],T[1]^2) ) and
 T[2]>=s2 and T[2]<=t2
then return true; fi;

if T[2]>=0 and 
(not HAP_SqrtStrictInequality(ABI*t2^2,T[3],ABI*T[2]^2) ) and
 T[1]>=s1 and T[1]<=t1
then return true; fi;


if T[1]<=0 and #=
(not  HAP_SqrtStrictInequality(s1^2,T[3],T[1]^2)  ) and
 T[2]>=s2 and T[2]<=t2
then return true; fi;


if T[2]<=0 and   #=
(not  HAP_SqrtStrictInequality(ABI*t2^2,T[3],ABI*T[2]^2)  ) and
 T[1]>=s1 and T[1]<=t1
then return true; fi;


if (T[1]-s1)^2 + ABI*(T[2]-s2)^2<=T[3]
then return true; fi;

if (T[1]-s1)^2 + ABI*(T[2]-t2)^2<=T[3]
then  return true; fi;

if (T[1]-t1)^2 + ABI*(T[2]-s2)^2<=T[3]
then  return true; fi;

if (T[1]-t2)^2 + ABI*(T[2]-t2)^2<=T[3]
then  return true; fi;

return false;
                                             #####
end);                                        #####
##################################################
##################################################


###################################################
###################################################
InstallGlobalFunction(IsQQUnimodularPair,     #####
function(arg)                                 #####
local OQ,x,BI,ABI,L,s1,t1,s2,t2,T,S;
#Inputs a unimodular pair x=[a,b] and returns true
#if the associated hemisphere touches the 
#fundamental Q-rectangle [0,t1]x[0,t2] in the complex 
#plane.


OQ:=arg[1];
BI:=OQ!.bianchiInteger;
ABI:=AbsInt(BI);
x:=arg[2];
L:=HAP_BianchiFundamentalRectangle(OQ);
s1:=0;t1:=L[2];; s2:=0; t2:=L[4];
#s1:=L[1];t1:=L[2];; s2:=L[3]; t2:=L[4];

T:=UnimodularPairCoordinates(OQ,x);

if T[1]>=0 and   #=
(not  HAP_SqrtStrictInequality(t1^2,T[3],T[1]^2) )
and T[2]>=0 and T[2]<=t2
then return true; fi;

if T[2]>=0 then    #=
if (not HAP_SqrtStrictInequality(ABI*t2^2,T[3],ABI*T[2]^2) )
and T[1]>=0 and T[1]<=t1
then return true; fi;
fi;

if T[1]<=0 then #=
if T[1]^2<=T[3]
and T[2]>=0 and T[2]<=t2
then return true; fi;
fi;

if T[2]<=0 then  #=
if ABI*T[2]^2<=T[3]
and T[1]>=0 and T[1]<=t1
then return true; fi;
fi;

if (T[1]-s1)^2 + ABI*(T[2]-t1)^2<=T[3]
then return true; fi;

if (T[1]-s1)^2 + ABI*(T[2]-t2)^2<=T[3]
then return true; fi;

if (T[1]-s2)^2 + ABI*(T[2]-t1)^2<=T[3]
then return true; fi;

if (T[1]-s2)^2 + ABI*(T[2]-t2)^2<=T[3]
then return true; fi;

return false;
                                             #####
end);                                        #####
##################################################
##################################################


################################################################
################################################################
InstallGlobalFunction(SimplicialComplexOfUnimodularPairs,  #####
function(OQ,L)                                             #####
local E,S,K,Lines,Points,bool,p,i,j,k,x;                  
#Inputs a list of unimodular pairs corresponding to disks 
#in the complex plane. The nerve of these disks is
#returned as an abstract simplicial complex.

Lines:=[];
K:=List(L,x->[Position(L,x)]);;
E:=[];
for i in [1..Length(L)] do
Lines[i]:=[];
for j in [i+1..Length(L)] do
if HAP_AreIntersectingUnimodularPairs(OQ,L[i],L[j])
then
Add(E,[i,j]);
fi;
od;od;

E:=SSortedList(E);

Append(K,E);

Points:=[];
for x in E do
i:=x[1]; j:=x[2];
for k in [j+1..Length(L)] do
if [i,k] in E and [j,k] in E then
  if not IsBound(Lines[i][j]) then
     Lines[i][j]:=UnimodularIntersectingLine(OQ,L[i],L[j]);
  fi;
  if not IsBound(Lines[i][k]) then
     Lines[i][k]:=UnimodularIntersectingLine(OQ,L[i],L[k]);
  fi;
  if not IsBound(Lines[j][k]) then
     Lines[j][k]:=UnimodularIntersectingLine(OQ,L[j],L[k]);
  fi;
  p:= HAP_Are3IntersectingUnimodularPairs(OQ,L[i],L[j],L[k],
                                  Lines[i][j],Lines[i][k],Lines[j][k]);
  if p[1] then
     Add(K,[i,j,k]);
     if not IsBound(p[2]) then
         p[2]:=UnimodularIntersectingPoint(OQ,
                  L[i],L[j],L[k],Lines[i][j],Lines[i][k]);
     fi;
     Add(Points,p[2]);
  fi;
fi;
od;od;


K:=MaximalSimplicesToSimplicialComplex(K);
K!.Lines:=Lines;
K!.Points:=Points;
K!.unimodularPairs:=L;
return K;                                      #####
end);                                          #####
####################################################
####################################################


#####################################################
#####################################################
InstallGlobalFunction(QQNeighbourhoodOfUnimodularPairs,
function(arg)                                   #####
local OQ,L,U,F,FF, Q, BI, x, p, u, S,St;      #####
#Inputs a list of unimodular pairs and returns a list 
#of all pairs/hemispheres in L that intersect
#the rectangle Q=[0,t1]x[0,t2]. 
OQ:=arg[1];
L:=arg[2];
F:=[];
Q:= Filtered(L,x->IsQQUnimodularPair(OQ,x));   
Q:=List(Q,x->UnimodularPairStandardForm(x));
BI:=OQ!.bianchiInteger;
if -BI mod 4 = 3 then
St:=[QuadraticNumber(1/2,1/2,BI),QuadraticNumber(-1/2,1/2,BI)];
else
St:=[QuadraticNumber(0,1,BI)];
fi;

for x in Q do
for S in St do
Add(F,x);
Add(F,[QuadraticNumberConjugate(x[1]),QuadraticNumberConjugate(x[2]),x[3]]);
Add(F,[x[1],-x[2],x[3]]);
Add(F,[QuadraticNumberConjugate(x[1]),-QuadraticNumberConjugate(x[2]),x[3]]);
##############
Add(F,[x[1],x[2]-x[1],x[3]]);
Add(F,[QuadraticNumberConjugate(x[1]),-QuadraticNumberConjugate(x[2])+QuadraticNumberConjugate(x[1]),x[3]]);
Add(F,[QuadraticNumberConjugate(x[1]),QuadraticNumberConjugate(x[2])-QuadraticNumberConjugate(x[1]),x[3]]);
Add(F,[x[1],-x[2]+x[1],x[3]]);
##############
Add(F,[x[1],x[2]-S*x[1],x[3]]);
Add(F,[QuadraticNumberConjugate(x[1]),-QuadraticNumberConjugate(x[2])-S*QuadraticNumberConjugate(x[1]),x[3]]);
Add(F,[x[1],x[2]-(1+S)*x[1],x[3]]);
Add(F,[QuadraticNumberConjugate(x[1]),-QuadraticNumberConjugate(x[2])+(1-S)*QuadraticNumberConjugate(x[1]),x[3]]);
##############
Add(F,[QuadraticNumberConjugate(x[1]),QuadraticNumberConjugate(x[2])+S*QuadraticNumberConjugate(x[1]),x[3]]);
Add(F,[x[1],-x[2]+S*x[1],x[3]]);
Add(F,[QuadraticNumberConjugate(x[1]),QuadraticNumberConjugate(x[2])+(-1+S)*QuadraticNumberConjugate(x[1]),x[3]]);
Add(F,[x[1],-x[2]+(1+S)*x[1],x[3]]);
od;
od;

F:=Classify(F,x->UnimodularPairStandardForm(x){[1,2]});;
F:=List(F,x->x[1]);
F:=List(F,x->x{[1,2,3]});

U:=Filtered(F,x->IsQQUnimodularPair(OQ,x));     
FF:=[];
for x in F do
for u in U do
if HAP_AreIntersectingUnimodularPairs(OQ,x,u) then
Add(FF,x); break; fi;
od;od;

return FF;                                       #####
end);                                           #####
#####################################################
#####################################################


#####################################################
#####################################################
InstallGlobalFunction(QNeighbourhoodOfUnimodularPairs,
function(arg)                                   #####
local OQ,L,U,F,FF, Q, BI, x, p, u, S,St;      #####
#Inputs a list of unimodular pairs and returns a list 
#of the 'images' of pairs in L that intersect the  
#fundamental rectangle [s1,t1]x[s2,t2]. 
OQ:=arg[1];
L:=arg[2];
F:=[];
Q:= Filtered(L,x->IsQQUnimodularPair(OQ,x));   
Q:=List(Q,x->UnimodularPairStandardForm(x));
BI:=OQ!.bianchiInteger;
if -BI mod 4 = 3 then
St:=[QuadraticNumber(1/2,1/2,BI),QuadraticNumber(-1/2,1/2,BI)];
else
St:=[QuadraticNumber(0,1,BI)];
fi;

for x in Q do
for S in St do
Add(F,x);
Add(F,[QuadraticNumberConjugate(x[1]),QuadraticNumberConjugate(x[2]),x[3]]);
Add(F,[x[1],-x[2],x[3]]);
Add(F,[QuadraticNumberConjugate(x[1]),-QuadraticNumberConjugate(x[2]),x[3]]);
##############
Add(F,[x[1],x[2]-x[1],x[3]]);
Add(F,[QuadraticNumberConjugate(x[1]),-QuadraticNumberConjugate(x[2])+QuadraticNumberConjugate(x[1]),x[3]]);
Add(F,[QuadraticNumberConjugate(x[1]),QuadraticNumberConjugate(x[2])-QuadraticNumberConjugate(x[1]),x[3]]);
Add(F,[x[1],-x[2]+x[1],x[3]]);
##############
Add(F,[x[1],x[2]-S*x[1],x[3]]);
Add(F,[QuadraticNumberConjugate(x[1]),-QuadraticNumberConjugate(x[2])-S*QuadraticNumberConjugate(x[1]),x[3]]);
Add(F,[x[1],x[2]-(1+S)*x[1],x[3]]);
Add(F,[QuadraticNumberConjugate(x[1]),-QuadraticNumberConjugate(x[2])+(1-S)*QuadraticNumberConjugate(x[1]),x[3]]);
##############
Add(F,[QuadraticNumberConjugate(x[1]),QuadraticNumberConjugate(x[2])+S*QuadraticNumberConjugate(x[1]),x[3]]);
Add(F,[x[1],-x[2]+S*x[1],x[3]]);
Add(F,[QuadraticNumberConjugate(x[1]),QuadraticNumberConjugate(x[2])+(-1+S)*QuadraticNumberConjugate(x[1]),x[3]]);
Add(F,[x[1],-x[2]+(1+S)*x[1],x[3]]);
od;
od;

F:=Classify(F,x->UnimodularPairStandardForm(x){[1,2]});;
F:=List(F,x->x[1]);
F:=List(F,x->x{[1,2,3]});
#return F;    

U:=Filtered(F,x->IsQUnimodularPair(OQ,x));

return U;                                       #####
end);                                           #####
#####################################################
#####################################################



#####################################################
#####################################################
InstallGlobalFunction(HAP_VertexHeights,        #####
function(OQ,LL)                                 #####
local L,Lines, Points, Heights, Cusps, K, 
 s, c, D,i,bool,S,ABI,a,b;

ABI:=AbsInt(OQ!.bianchiInteger);
L:=QQNeighbourhoodOfUnimodularPairs(OQ,LL);    
K:=SimplicialComplexOfUnimodularPairs(OQ,L);  
Lines:=K!.Lines;
Points:=K!.Points;
S:=K!.simplicesLst[3];

Heights:=[];
for i in [1..Length(S)]  do
if Points[i] = fail then Heights[i]:=-infinity; else
Heights[i]:= HAP_HeightOfPointOnSphere(OQ,Points[i],L[S[i][1]]);
fi;
od;
K!.Heights:=Heights;

Cusps:=[];
for i in [1..Length(Points)] do
if Heights[i]=0 then
bool:=true;
for s in L do 
c:=UnimodularPairCoordinates(OQ,s);
b:=Points[i][2]*(1/QuadraticNumber(0,1,ABI));
a:=(c[1]-Points[i][1])^2 + ABI*(c[2]-b)^2;
if a < c[3] then
bool:=false; break; fi;
od;
if bool then Add(Cusps, Points[i]); fi;
fi;
od;

Cusps:=SSortedList(Cusps);
K!.Cusps:=Cusps;

return [Cusps,K];                               #####
end);                                           #####
#####################################################
#####################################################


###################################################
###################################################
InstallGlobalFunction(UnimodularPairsReduced_NN, #####
function(OQ,LL)                                  #####
local L, Lines, Points, K, K3, E,EE, V, VV,
fn, fn2, v, w, h, dist, bool, BOOL, C,D,DD,i,F,S,M;

#L:=QQNeighbourhoodOfUnimodularPairs(OQ,LL);
if LL=[] then return [LL,[]]; fi;
L:=LL;

K:=SimplicialComplexOfUnimodularPairs(OQ,L);
Lines:=K!.Lines;
Points:=K!.Points;

K3:=K!.simplicesLst[3];

DD:=Filtered([1..K!.nrSimplices(2)], v-> IsQQUnimodularPair(OQ,L[K3[v][1]])
or IsQQUnimodularPair(OQ,L[K3[v][2]])
or IsQQUnimodularPair(OQ,L[K3[v][3]]) );
#Need to check that this is OK here

S:=K3{DD};
Points:=Points{DD};

V:=List([1..Length(S)],x->[S[x][1],S[x][2],S[x][3], Points[x] ] );

V:=Filtered(V,v->not v[4]=fail);

V:=List(V,v-> [v[1],v[2],v[3],v[4],HAP_HeightOfPointOnSphere(OQ,v[4],L[v[1]])]);
V:=Filtered(V,v->v[5]>=0);

VV:=[];
for v in V do
#h:=HAP_HeightOfPointOnSphere(OQ,v[4],L[v[1]]) ;
h:=v[5];
bool:=true; ;
#NEXT LINE IS TIME CONSUMING
#DD:=Filtered(K!.simplicesLst[3],x-> (v[1] in x and v[2] in x) or (v[1] in x and v[3] in x) or (v[3] in x and v[2] in x));  #?????????
DD:=Filtered(K!.simplicesLst[3],x-> (v[1] in x and v[2] in x) );

DD:=Concatenation(DD);
DD:=SSortedList(DD);
for i in DD do
w:=L[i];
#for w in L do
if  HAP_HeightOfPointOnSphere(OQ,v[4],w) - h > 0  then
bool:=false; break; fi;
od;
if bool then Add(VV,v); fi;
od;



V:=List(VV,x->[SortedList([x[1],x[2],x[3]]),x[4],x[5]]);
D:=K!.vertices; ##??
DD:=[];
for v in D do
F:=Filtered([1..Length(V)],i-> v in V[i][1]);
F:=List(F,i->[V[i][2],V[i][3]]);
F:=SSortedList(F);
if Length(F)>2 then Add(DD,v);  fi;
od;

L:=L{DD};
L:=QQNeighbourhoodOfUnimodularPairs(OQ,L); 
L:=SSortedList(L);


VV:=Filtered(VV,v->v[1] in DD and v[2] in DD and v[3] in DD);
M:=List(VV,v->Position(K!.simplicesLst[3],[v[1],v[2],v[3]]));
M:=Difference([1..Length(K!.Points)], M);
for i in M do
K!.Points[i]:=fail;
#K!.Heights:=-infinity;
od;


return [L,K];                                    #####
end);                                            #####
##################################################
##################################################

####################################################
####################################################
InstallGlobalFunction(UnimodularPairsReduced,
function(OQ,LL)
local K,L;

L:=QNeighbourhoodOfUnimodularPairs(OQ,LL);
if OQ!.bianchiInteger=-5 then
L:=Filtered(L,x->IsQUnimodularPair(OQ,x));
else
L:=Filtered(L,x->IsQQUnimodularPair(OQ,x));
fi;
K:=UnimodularPairsReduced_NN(OQ,L);
K:=QQNeighbourhoodOfUnimodularPairs(OQ,K[1]);
K:=UnimodularPairsReduced_NN(OQ,K);
return K;
end);
####################################################
####################################################

####################################################
####################################################
InstallGlobalFunction(CoverOfUnimodularPairs,  #####
function(arg)                                  #####
local OQ,N, NRMS, L, LL, K, Y, bool, BI, A,     
ABI, LST, UU, u, w, a, b, I, i, pos, HAPRECORD;

OQ:=arg[1];
N:=1;
BI:=OQ!.bianchiInteger;
ABI:=AbsInt(BI);
###############################

HAPRECORD:=[
[ 1, 1 ],[ 2, 1 ],[ 3, 1 ],[ 5, 20 ],[ 6, 24 ],[ 7, 1 ],[ 10, 40 ],[ 11, 1 ],[ 13, 52 ],[ 14, 56 ],[ 15, 15 ],[ 17, 68 ],[ 19, 5 ],[ 21, 84 ],[ 22, 88 ],[ 23, 18 ],[ 26, 121 ],[ 29, 169 ],[ 30, 169 ],[ 31, 20 ],[ 33, 528 ],[ 34, 225 ],[ 35, 35 ],[ 37, 289 ],[ 38, 4 ],[ 39, 39 ],[ 41, 361 ],[ 42, 672 ],[ 43, 9 ],
#[ 46, 25 ],
[ 47, 36 ],[ 51, 51 ],
#[ 53, 25 ],
[ 55, 56 ],
#[ 57, 49 ],
[ 59, 36 ],[ 67, 23 ],[ 71, 64 ],[ 79, 64 ],[ 83, 36 ],[ 87, 87 ],[ 91, 91 ],[ 95, 95 ],[ 103, 64 ],[ 107, 64 ],[ 111, 111 ],[ 115, 115 ],
#[ 119, 25 ],
[ 123, 123 ],[ 127, 121 ],[ 131, 100 ],[ 139, 100 ],[ 143, 576 ],[ 151, 110 ],[ 155, 155 ],[159,240],[ 163, 53 ],[167,384],[ 179, 144], [199,196]
];;
                 #HAPRECORD contains some precomuted [d,r] for which Swan's
                 #criterion is satisfied for the given norm r.
#HAPRECORD:=[];  #Uncomment this line to calculate HAPRECORD entries from 
                 #scratch.  

if Length(arg)=2 then pos:=arg[2];
else

pos:=Position(List(HAPRECORD,x->x[1]),ABI);
if pos =fail or pos=0 then pos:=infinity;
if Length(arg)=1 then
Print("Try \n  P:=BianchiPolyhedron(OQ,N);\nfor some guessed positive integer value of N and then try\n  SwanBianchiCriterion(P);\nto test if the value of N was large enough. If the test returns false then you'll need to try a larger value of N.\n\nA successful value of N can be stored as a pair [d,N] in the list HAPRECORD which can be edited manually in the file hap/lib/Congruence/bianchi.gi .\n\n");
fi;
else pos:=HAPRECORD[pos][2]; fi;

fi;
###############################
NRMS:=QuadraticIntegersByNorm(OQ,1000);  #We'll never need anything higher!!!
NRMS:=List(NRMS,x->HAPNorm(OQ,x));
NRMS:=SSortedList(NRMS);
bool:=true;

L:=[];

while bool do
N:=N+1; 
if Length(arg)>1 then
#Print("Adding hemispheres of squared radius ",1/NRMS[N],"\n");
fi;

L:=UnimodularPairs(OQ,NRMS[N],true,L);
if pos=infinity then
K:=SimplicialComplexOfUnimodularPairs(OQ,
                 QNeighbourhoodOfUnimodularPairs(OQ,L));
Y:=RegularCWComplex(K);
CocriticalCellsOfRegularCWComplex(Y,1);
if Homology(Y,0)=[0] and Homology(Y,1)=[] then break; fi;
fi;
if NRMS[N]>=pos then break; fi;
od;

K:=UnimodularPairsReduced(OQ,L);

L:=K[1]; K:=K[2];

###########CHECK SWAN'S CONDITION REALLY IS MET
if pos=infinity then 
A:=BianchiPolyhedron(OQ!.bianchiInteger,L);
A:=SwanBianchiCriterion_alt(A);
if not Length(A)=0 then 
return CoverOfUnimodularPairs(OQ,N+1,true);
fi;
fi;
##################################################

L:=Filtered(L,x->IsQQUnimodularPair(OQ,x));   
L:=QQNeighbourhoodOfUnimodularPairs(OQ,L);
L:=Filtered(L,x->IsQQUnimodularPair(OQ,x));
return [L,K];                                    #####
end);                                        #####
##################################################
##################################################

##################################################
##################################################
InstallGlobalFunction(IsUnimodularCollection,
function(OQ,L)
local K,Y;
K:=SimplicialComplexOfUnimodularPairs(OQ,L);
Y:=RegularCWComplex(K);
CocriticalCellsOfRegularCWComplex(Y,1);
if not Homology(Y,0)=[0] then return false; fi;
if not Homology(Y,1)=[] then return false; fi;
return true; 
end);
##################################################
###############################################

##################################################
##################################################
InstallGlobalFunction(SwanBianchiCriterion,
function(P)
local mr, mv, R, L,LL, i, u, Points;

mr:=P!.minRadius;
mv:=P!.minVertexHeight;
if IsInt(1/mv) then R:=1/mv;
else R:= 1+Int(1/mv);
fi;

L:=UnimodularPairs(P!.ring,R);;
LL:=UnimodularPairs(P!.ring,mr);;
L:=Difference(L,LL);;
L:=Filtered(L,x->IsQQUnimodularPair(P!.ring,x));;
L:=Filtered(L,x->not HAP_IsRedundantUnimodularPair(P!.ring,P!.unimodularPairs,x));

Points:=Filtered([1..Length(P!.points)],i->P!.visibleHeights[i]>0);

for u in L do
for i in Points do
if P!.visibleHeights[i] < HAP_HeightOfPointOnSphere(P!.ring,P!.points[i],u) then return false; fi;
od;
od;

return true;

end);
##################################################
##################################################


##################################################
##################################################
InstallGlobalFunction(SwanBianchiCriterion_alt,
function(arg)
local  PP,P,F,A,K,C,CC,D,H, nrmL, L, S, S1, v,b,x,i,h,s1,s2,t1,t2;

P:=arg[1];

if not IsUnimodularCollection(P!.ring,P!.unimodularPairs) then 
return [false]; fi;

H:=P!.heights;
H:=Filtered(H,x->x>0);
H:=1/Minimum(H);
if IsHapQuadraticNumber(H) then H:=H!.rational; fi;

L:=P!.unimodularPairs;
L:=List(L,x->HAPNorm(P!.ring,x[1]));;
nrmL:=Maximum(L);

if nrmL>=H then return []; fi;

L:=UnimodularPairs(P!.ring,nrmL);

if Length(arg)=2 then 
   H:=UnimodularPairs(P!.ring,arg[2]);
   H:=QNeighbourhoodOfUnimodularPairs(P!.ring,H);
else
   H:=UnimodularPairs(P!.ring,H);
   H:=Filtered(H,h->not h in L);
   H:=QNeighbourhoodOfUnimodularPairs(P!.ring,H);
fi;


S:=QuadraticNumber(0,1,-P!.bianchiInteger);
S1:=1/S;

L:=HAP_BianchiFundamentalRectangle(P!.ring);
s1:=0; t1:=L[2]; s2:=0; t2:=L[4];
#s1:=L[1];t1:=L[2];; s2:=L[3]; t2:=L[4];


D:=[];
for i in [1..Length(P!.points)] do
if not P!.points[i]=fail then
if P!.heights[i]>0 
and P!.points[i][1]>=s1 and P!.points[i][1]<=t1 and P!.points[i][2]*S1>=s2
and P!.points[i][2]*S1<=t2
then Add(D,i);
fi;
fi;
od;

C:=[]; CC:=[];
for i in D do
for h  in H do
if HAP_HeightOfPointOnSphere(P!.ring,P!.points[i],h) > P!.heights[i]
then 
Add(C,h);
Add(CC,i);
fi;
od;
od;
CC:=SSortedList(CC);

K:=P!.simplicialComplex;
A:=[];
for v in K!.vertices do
F:=Filtered(K!.simplicesLst[3],x-> v in x);
F:=Concatenation(F);
F:=SSortedList(F);
if IsSubset(CC,F) then Add(A,K!.unimodularPairs[v]); fi;
od;
if Length(A)>0 then return A; fi;

C:=Filtered(C,x->x[3]>nrmL);
PP:=BianchiPolyhedron(P!.bianchiInteger,Concatenation(P!.unimodularPairs,C));
return Filtered(PP!.unimodularPairs,x->not x in P!.unimodularPairs);

end);
##################################################
##################################################

#####################################################
#####################################################
InstallGlobalFunction(BianchiPolyhedron,
function(arg)
local d,K, P,C, OQ, F, H, L,M,i,x,N,s;

d:=arg[1];
if d>=0 then
Print("Input must be a negative square free integer.\n");
return fail;
fi;

F:=Factors(AbsInt(d));

if not Length(F)=Length(SSortedList(F)) then
Print("Input is not a square free integer.\n");
return fail;
fi;

P:=Objectify(HapBianchiPolyhedron,rec());
P!.bianchiInteger:=d;
OQ:=RingOfIntegers(QuadraticNumberField(d));
P!.ring:=OQ;

if AbsInt(d) in [1,3] then
Print("Warning: the polyhedron is a union of fundamental domains for PSL( ",OQ, ").\n");
fi; 

if Length(arg)=1 then
    K:=CoverOfUnimodularPairs(OQ); P!.unimodularPairs:=K[1]; K:=K[2];
else

    if IsInt(arg[2]) then
        K:=CoverOfUnimodularPairs(OQ,arg[2]); P!.unimodularPairs:=K[1]; K:=K[2];
    else
        K:=UnimodularPairsReduced(OQ,arg[2]); P!.unimodularPairs:=K[1]; K:=K[2];
    fi;
fi;

N:=4;
while Length(P!.unimodularPairs)=0 do 
#Print("Try again with a larger choice of limit on the norms or a larger list of unimodular pairs\n");
N:=N+1;     #Bit of a botch job!
K:=CoverOfUnimodularPairs(OQ,N); P!.unimodularPairs:=K[1]; K:=K[2];
od;


K:=HAP_VertexHeights(OQ,P!.unimodularPairs);

P!.unimodularPairs:=K[2]!.unimodularPairs;
P!.simplicialComplex:=K[2];
P!.vertices:=K[2]!.Points;
P!.heights:=K[2]!.Heights;
C:=[];
for x in K[1] do
Add(C,x);
Add(C,[-x[1],x[2]]);
Add(C,[x[1],-x[2]]);
Add(C,[-x[1],-x[2]]);
od;
C:=SSortedList(C);
P!.cusps:=C;




L:=K[2]!.unimodularPairs;
H:=P!.heights;
############
for i in [1..Length(H)] do
if H[i]>0 then
for s in L do 
if HAP_HeightOfPointOnSphere(OQ,K[2]!.Points[i],s)!.rational>H[i] then
H[i]:=-infinity;
fi;
od;
fi;
od;
############
P!.visibleHeights:=1*H;
H:=Filtered(H,x->x>0);
if Length(H)>0 then
H:=Minimum(H);
else
H:=0;
fi;

L:=P!.unimodularPairs;
L:=List(L,x->HAPNorm(P!.ring,x[1]));;
L:=1/Maximum(L);

P!.minVertexHeight:=H;
P!.minRadius:=L;
P!.points:=K[2]!.Points;

P!.unimodularPairs:=QNeighbourhoodOfUnimodularPairs(P!.ring,P!.unimodularPairs);
return P;

end);
####################################################
####################################################