######### PACKAGE FOR FINDING DIVISIBLE HYPERBOLIC QUADRILATERALS 

with(linalg):
with(numtheory):

## GEOMETRIC FUNCTIONS
# equation solving functions

FilterComplexSolutions := proc(SL,eps)
   local LSL,sol,csol,eq,var,val,valid;
   LSL := [];
   for sol in SL do
       valid := true;
       csol := {};
       for eq in sol do
           var := lhs(eq);
           val := rhs(eq);  
           if evalc(abs(Im(val))) > eps 
              then valid := false;
              else val :=  evalc(Re(val));
              fi;
           csol := {op(csol),var=val};
           od;
       if valid then LSL :=[op(LSL),csol]; fi;
       od;
   RETURN(LSL);
   end:         
       

# basic hyperbolic geometry constructs

EuctoHyp := proc(r)
   RETURN(ln((1+r)/(1-r)));
   end:

HyptoEuc := proc(h)
   RETURN(tanh(h/2));
   end:

hdist := proc(P,Q)
   local z0,z1,bz0,ab;
   z0 :=P[1]+P[2]*I;
   z1 := Q[1]+Q[2]*I;;
   bz0 := evalc(conjugate(z0));
   ab :=  normal(evalc(abs((z1-z0)/(1-bz0*z1))));
   RETURN(evalf(EuctoHyp(ab)));
   end:

Zhline := proc(x0,y0,x1,y1)
           local eq0,eq1,a,r,Sol;
           eq0 := (x0-a)^2+(y0-b)^2=a^2+b^2-1;
           eq1 := (x1-a)^2+(y1-b)^2=a^2+b^2-1;
           Sol := solve({eq0,eq1},{a,b});
           RETURN(simplify(subs(Sol,[a,b])));
           end:
 
hlinecirc := proc(x0,y0,x1,y1,t)
         local a,b,r,Z;
         Z :=Zhline(x0,y0,x1,y1);
         a :=Z[1]; b := Z[2];
         r := sqrt(a^2+b^2-1);
         RETURN([a+r*cos(t),b+r*sin(t),t=0..2*Pi]);
         end:
 
IntPoint := proc(a1,b1,a2,b2)
    local eq1,eq2,Sol,P,r1,r2,r3,R;
    r1 := sqrt(a1^2+b1^2-1);
    r2 := sqrt(a2^2+b2^2-1);
    r3 := sqrt((a1-a2)^2+(b1-b2)^2);
    if ((r1 >= r2+r3) or
        (r2 >= r1+r3) or
        (r3 >= r1+r2))  
        then RETURN(FAIL); fi;
    eq1 := 2*a1*x+2*b1*y = 1+x^2+y^2;
    eq2 := 2*a2*x+2*b2*y = 1+x^2+y^2;
    eq2 := lhs(eq1)-lhs(eq2)=0; 
    Sol := [fsolve({eq1,eq2},{x,y},{x=-100..100,y=-1000..100})];
    if nops(Sol) = 0 then RETURN(FAIL); fi;
    P := subs(Sol[1],[x,y]);
    R := P[1]^2+P[2]^2;
    if R > 1 then P := [P[1]/R,P[2]/R]; fi;
    RETURN(P);
    end:

PointandAngle:= proc(a1,b1,a2,b2)

     local P,r1,r2,r3,C;
     P := IntPoint(a1,b1,a2,b2);
     if P = FAIL then RETURN(FAIL); fi;
     r1 := sqrt((a1-P[1])^2+(b1-P[2])^2);
     r2 := sqrt((a2-P[1])^2+(b2-P[2])^2);
     r3 := sqrt((a1-a2)^2+(b1-b2)^2);
     C  := evalf(arccos((r3^2-r1^2-r2^2)/(2*r1*r2)));
     RETURN([P,C]);
     end:


CentreAngleEqn := proc(A,B,a,b,C)
   local eq,R,r;
   eq := (A-a)^2+(B-b)^2 = R^2+r^2+2*R*r*cos(C); 
   eq := simplify(subs(R^2=A^2+B^2-1,r^2=a^2+b^2-1,eq));
   if cos(C) = 0 then RETURN(eq); fi;  
   eq := R*r=simplify(r*rhs(isolate(eq,R)));
   eq := lhs(eq)^2 = rhs(eq)^2;
   eq := subs(R^2=A^2+B^2-1,r^2=a^2+b^2-1,eq);
   eq := expand(eq);
   RETURN(lhs(eq)-rhs(eq)=0);
   end:

Transdist :=  proc(A)
   local d,t,lam;
   d:=det(A);
   t := trace(A);
   t := abs(t/sqrt(d));
   if (t < 2) 
        then ERROR(`the element is not hyperbolic`);
        else RETURN(2*arccosh(t/2));
        fi;           
     end:


#schwartz triangle functions

SchwartzCenter:= proc(a,b,c)
    local ca,sa,cb,sb,cc,sc,H1,H2,AB,x,y,z,bz;
    ca := cos(a); sa :=sin(a);
    cb := cos(b); sb :=sin(b);
    cc := cos(c); sc :=sin(c); 
    H1 :=cb^2*x^2-sb^2*y^2-cb^2; 
    H1 := evalf(H1); 
    H2 := cc^2*x^2-sc^2*y^2-cc^2;
    H2 := subs({x=ca*x+sa*y,y=-sa*x+ca*y},H2); 
    H2 := evalf(H2);          
    AB:= fsolve({evalf(H1=0),evalf(H2=0)},{x,y}, x=0..1000,y=0..1000);
    RETURN(subs(AB, x+y*I)); 
    end:

SchwartzT := proc(a1,a2,a3)
     local l,m,n,z,a,b,x0,u,s,t,x,y,P1,P2,P3; 
     z := SchwartzCenter(a1,a2,a3);
     a := evalc(Re(z));
     b := evalc(Im(z));
     P1 := [0,0];
     P2 := [fsolve(x^2-2*a*x+1,x,-1..1),0];
     P3 := [fsolve({2*a*x+2*b*y=1+x^2+y^2,
                    sin(a1)*x-cos(a1)*y=0},
                    {x,y},x=-1..1,y=0..1)];
     P3 := subs(P3[1],[x,y]);
     RETURN([P1,P2,P3]);
     end:

SchwartzSides := proc(a1,a2,a3)
   local pnts,P1,P2,P3;  
   pnts := SchwartzT(a1,a2,a3);
   P1 := pnts[1];
   P2 := pnts[2];
   P3 := pnts[3];
   RETURN([hdist(P1,P2),hdist(P2,P3),hdist(P3,P1)]);
   end:

DrawSchwartzT:= proc(a1,a2,a3)
  local pnts,P1,P2,P3,t,line1,line2,circ0,circ1;  
  pnts := SchwartzT(a1,a2,a3);
  P2 := pnts[2];
  P3 := pnts[3];
  line1 := [t,0,t=-1..1];
  line2 := [t*cos(a1),t*sin(a1),t=-1..1];
  circ0 := [cos(t),sin(t),t=0..2*Pi]; 
  circ1 := hlinecirc(P2[1],P2[2],P3[1],P3[2],t);
  plot({circ0,line1,line2,circ1},view=[-1.2..1.2,-1.2..1.2],
        color=black, axes=none,scaling=constrained);
  end:
  


# Lmn matrix triples


FlmnTriple:= proc(lmn)
    local l,m,n,cl,sl,z0,bz0,Pm,Qm,Rm,Pmc,Qmc,Rmc,Am,Bm,Cm;
    l := lmn[1]; m :=lmn[2]; n := lmn[3]; 
    cl := cos(Pi/l); sl :=sin(Pi/l);
    z0 := SchwartzCenter(Pi/l,Pi/m,Pi/n);
    bz0 := conjugate(z0);
    Pm := matrix(2,2,evalf([cl+I*sl,0,0,cl-I*sl]));
    Qm := matrix(2,2,[1,0,0,1]);
    Rm := matrix(2,2,[z0,-1,1,-bz0]);
    Pm := evalm(sqrt(det(Pm))^(-1)*Pm);
    Qm := evalm(sqrt(det(Qm))^(-1)*Qm);
    Rm := evalm(sqrt(det(Rm))^(-1)*Rm);
    Pmc := htranspose(transpose(Pm));
    Qmc := htranspose(transpose(Qm));
    Rmc := htranspose(transpose(Rm));
    Am :=evalm(Pm&*Qmc);
    Bm :=evalm(Qm&*Rmc);
    Cm :=evalm(Rm&*Pmc);
    RETURN([Am,Bm,Cm]);
    end:

IsPSL2Identity:=  proc(A,eps)
    local B,Id,err,i,j;
    B := evalm(A);
    if evalf(abs(B[1,1])) <> 0 then B := evalm((1/B[1,1])*B); fi;
    Id := matrix(2,2,[1,0,0,1]);
    err := evalm(B-Id);
    for i from 1 to 2 do
    for j from 1 to 2 do
        if evalf(abs(err[i,j])) > eps then RETURN(false); fi:
        od;
        od;
    RETURN(true);
    end:


# hyperbolic quadilateral functions

MakeQuad := proc(a1,a2,a3,a4,len)
    local a,b,H,L,A1,A2,B1,B2,A3,B3,AB,eq1,eq2,x0,
          P0,P1,P2,P3,P4,eps,x,y,C1,C2,sol;
    H := cos(a1)^2*a^2-sin(a1)^2*b^2 = cos(a1)^2; 
    x0 := HyptoEuc(len/2);
    L := 2*a*x+2*b*y=1+x^2+y^2;
    L := subs(x=-x0,y=0,L);
    AB := fsolve({H,L},{a,b},a=-1000..1000,b=0..1000);
    A1 := subs(AB,a);  B1 := subs(AB,b);
    H := cos(a2)^2*a^2-sin(a2)^2*b^2 = cos(a2)^2; 
    L := 2*a*x+2*b*y=1+x^2+y^2;
    L := subs(x=x0,y=0,L);
    AB := fsolve({H,L},{a,b},a=0..1000,b=0..1000);
    A2 := subs(AB,a);  B2 := subs(AB,b);
    P0 := IntPoint(A1,B1,A2,B2);
    if (not ( P0 = FAIL )) then RETURN(FAIL); fi;
    eq1 := CentreAngleEqn(A1,B1,a,b,a4);
    eq2 := CentreAngleEqn(A2,B2,a,b,a3);
    AB := FilterComplexSolutions([solve(evalf({eq1,eq2}),{a,b})],_gGeoeps);
    A3 := FAIL;
    eps := max((0.1)*min(abs(evalf(Pi-2*a3)),abs(evalf(Pi-2*a4))),10^(2-Digits/2));
    for sol in AB do
        x := subs(sol,a); 
        y := subs(sol,b); 
        C1 := PointandAngle(A1,B1,x,y);
        C2 := PointandAngle(A2,B2,x,y); 
        if  (not (C1=FAIL or C2=FAIL))
            then if ((abs(evalf(C1[2]-a4))+abs(evalf(C2[2]-a3)) < eps)
                     and (C1[1][2] >0) 
                     and (C2[1][2] >0)) 
                     then A3 := x; B3 := y;
                     fi;
            fi;
        od;
    if A3 = FAIL then RETURN(FAIL); fi;
    P1 := [-x0,0];
    P2 := [x0,0];
    P3 := IntPoint(A2,B2,A3,B3);
    P4 := IntPoint(A1,B1,A3,B3);
    if ((P3 = FAIL) or (P4 = FAIL)) then RETURN(FAIL);fi;
    RETURN([P1,P2,P3,P4]);
    end: 

QuadLengths := proc(a1,a2,a3,a4,len)
  local Q,P1,P2,P3,P4;
  Q :=MakeQuad(a1,a2,a3,a4,len);
  if Q = FAIL then RETURN(FAIL); fi; 
  P1 := Q[1]; P2 :=Q[2]; P3 :=Q[3]; P4 :=Q[4];
  RETURN([hdist(P1,P2),hdist(P2,P3),hdist(P3,P4),hdist(P4,P1)]);
  end:

 DrawQuad := proc(a1,a2,a3,a4,len)
  local Q,P1,P2,P3,P4,line,circ0,circ1,circ2,circ3,t;
  Q :=MakeQuad(a1,a2,a3,a4,len);
  P1 := Q[1]; P2 :=Q[2]; P3 :=Q[3]; P4 :=Q[4];
  if Q = FAIL then RETURN(FAIL); fi; 
  line := [t,0,t=-1..1];
  circ0 := [cos(t),sin(t),t=0..2*Pi]; 
  circ1 := hlinecirc(P2[1],P2[2],P3[1],P3[2],t);
  circ2 := hlinecirc(P3[1],P3[2],P4[1],P4[2],t);
  circ3 := hlinecirc(P4[1],P4[2],P1[1],P1[2],t);
  plot({circ0,line,circ1,circ2,circ3},view=[-1.2..1.2,-1.2..1.2],
        color=black, axes=none,scaling=constrained);
  end:
 

QuadBase := proc(r,s)
     local l,m,n,z,a,b,x0,u; 
     if (r+s = Pi) then RETURN(0); fi;   
     z := SchwartzCenter(r,s,0);
     a := evalc(Re(z));
     b := evalc(Im(z));
     x0 := fsolve(2*a*x = 1+x^2,x,-1..1);
     u := EuctoHyp(x0);
     RETURN(evalf(u));
     end:


## COMBINATORIAL FUNCTIONS

# List functions 

numequal := proc(L,x)
   local count,y;
   count := 0;
   for y in L do if x=y then count := count+1; fi; od;
   RETURN(count);
   end: 

indicesequal := proc(L,x)
    local M,i;
    M := [];
    for i from 1 to nops(L) do if x=L[i] then M := [op(M),i]; fi; od;
    RETURN(M);
    end:

indicessuperdivide := proc(L,x)
    local M,i;
    M := [];
    for i from 1 to nops(L) do if (L[i] mod x) = 0 then M := [op(M),i]; fi; od;
    RETURN(M);
    end:


allint := proc(L)
   local s;
   for s in L do
       if (not type(s,integer)) then RETURN(false) fi;
       od;
   RETURN(true);
   end:
 
ListLT :=  proc(a,b)
 local m,i;
 m := min(nops(a),nops(b));
 for i from 1 to m do
     if a[i] < b[i] then RETURN(true);  fi;
     if a[i] > b[i] then RETURN(false); fi;
     od;
 if nops(a) < nops(b) 
    then RETURN(true); 
    else RETURN(false);
    fi;
 end:
 
 CyclicNormalForm :=  proc(L)
    local LL,n,i;
    n := nops(L);
    LL :=[seq([op(L[i..n]),op(L[1..i-1])],i=1..n)];
    LL := sort(LL,ListLT);
    RETURN(LL[1]);
    end:



# tiling  functions

lammunu := proc(lmn);
  RETURN([seq(iquo(lmn[i],2),i=1...nops(lmn))]);
  end:

lmnDivisors := proc(lmn)
  local i,dset;
  dset := {};
  for i from 1 to nops(lmn) do
     dset := {op(dset),op(divisors(lmn[i]))};
     od;
  dset := dset minus {1};
  RETURN(sort([op(dset)]));
  end: 

tilemu := proc(angs);
   RETURN(nops(angs)-2-sum(1/'angs[i]','i'=1..nops(angs)));
   end: 

# quadrilateral and subdivided construction functions 

MakeQuadList := proc(AL)
   local i,j,k,l,n,IL;
   n := nops(AL);
   IL := {seq(seq(seq(seq([AL[i],AL[j],AL[k],AL[l]],
          i=1..n),j=1..n),k=1..n),l=1..n)};
   IL := [op(IL minus {[2,2,2,2]})];
   IL := {seq(CyclicNormalForm(IL[i]),i=1..nops(IL))};
   IL := [op(IL)];
   IL := sort(IL,ListLT);
   RETURN(IL);
   end:

Ksize := proc(Q,lmn);
   RETURN(tilemu(Q)/tilemu(lmn));
   end: 
 
Ktest := proc(Q,lmn,mK)
    local K;
    K := tilemu(Q)/tilemu(lmn); 
    if (K >= mK and type(K,integer)) 
       then RETURN(true) 
       else RETURN(false);
       fi;
    end: 

QuadMult := proc(Q,lmn,assn)
   RETURN([seq(lmn[assn[i]]/Q[i],i=1..nops(Q))]);
   end:
 
NumHubsAssn := proc(Q,lmn,assn)
  local i,K,ch,chi,hubs,mult;
  mult := QuadMult(Q,lmn,assn); 
  hubs := []; 
  K := Ksize(Q,lmn); 
  for i from 1 to nops(lmn) do
      chi := indicesequal(assn,i); 
      ch := sum('mult[chi[i]]','i'=1..nops(chi));
      hubs := [op(hubs),(K-ch)/lmn[i]]; 
      od;
  RETURN(hubs);
  end:


KandHubTest := proc(Q,lmn,assn,mK)

   local K,hubs;

   if (not Ktest(Q,lmn,mK)) then RETURN(false); fi; 

   hubs := NumHubsAssn(Q,lmn,assn);
   if (not allint(hubs)) 
      then RETURN(false);
      else RETURN(true);

      fi;
   end: 


MakeSQuad := proc(Q,lmn,assn)
   local K,mult,hubs;
   K := Ksize(Q,lmn); 
   mult := QuadMult(Q,lmn,assn);
   hubs := NumHubsAssn(Q,lmn,assn);
   RETURN(table([ (_ID) =[Q,assn],(_Q)=Q, (_T)=lmn, (_assn)=ha,
                 (_K)=K, (_mult)=mult,(_hubs) = hubs]));
    
   end:
 
 
SimpleHubAssign := proc(Q,lmn,mK)
   local i,j,k,l,R1,R2,R3,R4,assn,assn2;	
   R1 := indicesequal(lmn,Q[1]);
   R2 := indicesequal(lmn,Q[2]);
   R3 := indicesequal(lmn,Q[3]);
   R4 := indicesequal(lmn,Q[4]);
   assn := [seq(seq(seq(seq([R1[i],R2[j],R3[k],R4[l]],
           i=1..nops(R1)),j=1..nops(R2)),k=1..nops(R3)),l=1..nops(R4))];
   assn2 := [];
   for k from 1 to nops(assn) do
      if KandHubTest(Q,lmn,assn[k],mK) 
         then assn2 := [op(assn2),assn[k]];
         fi;      
      od;
   RETURN(assn2);
   end:


HubAssign := proc(Q,lmn,mK)
   local i,j,k,l,R1,R2,R3,R4,assn,assn2;	
   R1 := indicessuperdivide(lmn,Q[1]);
   R2 := indicessuperdivide(lmn,Q[2]);
   R3 := indicessuperdivide(lmn,Q[3]);
   R4 := indicessuperdivide(lmn,Q[4]);
   assn := [seq(seq(seq(seq([R1[i],R2[j],R3[k],R4[l]],
           i=1..nops(R1)),j=1..nops(R2)),k=1..nops(R3)),l=1..nops(R4))];
   assn2 := [];
   for k from 1 to nops(assn) do
      if KandHubTest(Q,lmn,assn[k],mK) 
         then assn2 := [op(assn2),assn[k]];
         fi;      
      od;
   RETURN(assn2);
   end:



# Segment functions 

 SegInit := proc(lmn,st,fn)
    local seg;
    seg := array;
    seg[_verts] := [st,fn];
    seg[_O] := 1;
    seg[_hlen] := _gTSides[st,fn];
    RETURN(seg);
    end:
 
 SegNumVerts := proc(seg);
    RETURN(nops(seg[_verts]));
    end:
 
 SegNumEdges := proc(seg);
    RETURN(SegNumVerts(seg)-1);
    end:
 
 SegVert := proc(seg,n);
    RETURN(seg[_verts][n]);
    end:
 
 SegVerts := proc(seg);
    RETURN(seg[_verts]);
    end:
 
 SegStart := proc(seg);
    RETURN(SegVert(seg,1));
    end:
 
 SegEnd := proc(seg);
    RETURN(SegVert(seg,SegNumVerts(seg)));
    end:
 
 SegPenVert := proc(seg);
    RETURN(SegVert(seg,SegNumVerts(seg)-1));
    end:

 SegEndOrient := proc(seg)
   local lseg,P,Q,R,exp,nP,sgn;

   P := SegPenVert(seg);

   Q := SegEnd(seg);
   R := (3+Q-P) mod 3; 

   if R <>  1 then R := -1; fi;
   RETURN(R);

   end:

 
 SegHlen := proc(seg)
  RETURN(seg[_hlen]);
  end: 
 
 SegExtend := proc(seg,lmn)
    local lseg,P,Q,R,exp,nP,sgn;
    lseg := copy(seg);
    P := SegPenVert(seg);
    Q := SegEnd(seg);
    R := {1,2,3} minus {P,Q};
    R := R[1];
    exp := lmn[Q];
    if type(exp,odd) then nP := R; sgn := 1; else nP := P; sgn := -1;fi;
    lseg[_verts] := [op(SegVerts(seg)),nP];
    lseg[_O] := lseg[_O]*sgn;
    lseg[_hlen]:= lseg[_hlen]+_gTSides[Q,nP];
    RETURN(lseg);
    end:

SegHubNums := proc(seg)
   local Iverts;
   Iverts := SegVerts(seg);
   Iverts := Iverts[2..(nops(Iverts)-1)];
   RETURN([seq(numequal(Iverts,i),i=1..3)]); 
   end:


SegMaxlen := proc(SQuad);
   RETURN(SQuad[_K]+2);
   end:

SegMinHlen := proc(i,j);
   RETURN(_gQBase[i,j]);
   end:

SegFind := proc(SL,h,eps)
   local seg;
   for seg in SL do
     if (abs(h-SegHlen(seg)) < eps) then RETURN(seg) fi;
     od;
   RETURN(FAIL);
   end:


# side functions


  
MakeSides := proc(lmn,pat,max)
    local a,b,nx,pv,cpat,apat,cpats,apats,i;
    a := pat[1];
    b := pat[2];
    nx := (a mod 3) +1;
    pv := ((a-2) mod 3) +1;  
    cpat := SegInit(lmn,a,nx);
    apat := SegInit(lmn,a,pv);
    cpats := [];
    apats := [];
    for i from 1 to max do
        if SegEnd(cpat)= b then cpats := [op(cpats),cpat]; fi;
        if SegEnd(apat)= b then apats := [op(apats),apat]; fi;
        cpat := SegExtend(cpat,lmn);
        apat := SegExtend(apat,lmn); 
        od;
    RETURN([cpats,apats]); 
    end:
 
MakeAllSides := proc(lmn,max)
    local allsides,i,j,pats;
    allsides := table();
    for i from 1 to nops(lmn) do
    for j from 1 to nops(lmn) do
        pats := MakeSides(lmn,[i,j],max);
        allsides[i,j,1]  := pats[1];
        allsides[i,j,-1] := pats[2];
        od;
        od;
    RETURN(allsides); 
    end:
 
FilterSideList := proc(allsides,SQuad)
    local SA,i,j,ii,jj,iii,jjj,n,assn,mx,mn,s,ne,le,Q;

    assn := SQuad[_assn];  
    Q := SQuad[_Q];

    mx := SegMaxlen(SQuad);
    SA := array;   
    n := nops(assn);

    for i from 1 to n do

      j := (i mod n) +1;
      ii := assn[i];
      jj := assn[j];

      iii := Q[i];
      jjj := Q[j];
      SA[i,1] := [];
      SA[i,-1] := [];
      mn  := SegMinHlen(iii,jjj);     
      for s in allsides[ii,jj,1] do
          ne := SegNumEdges(s);
          le := SegHlen(s);
          if ((le >= mn) and (ne <= mx)) 
             then SA[i,1] := [op(SA[i,1]),s]; 
             fi;
          od; 
      for s in allsides[ii,jj,-1] do
          ne := SegNumEdges(s);
          le := SegHlen(s);
          if ((le >= mn) and (ne <= mx)) 
             then SA[i,-1] := [op(SA[i,-1]),s]; 
             fi;
          od; 
      od;
    RETURN(SA);
    end:


## SETUP FUNCTIONS

SetTSides := proc(lmn)
  local l,m,n,T,TS; 
  TS := table(symmetric);
  T :=  SchwartzSides(Pi/lmn[1],Pi/lmn[2],Pi/lmn[3]);
  TS[1,2] := T[1];
  TS[1,3] := T[3];
  TS[2,3] := T[2];
  RETURN(TS);
  end:

SetQuadBases := proc(lmn)
  local i,j,QB,lmnD; 
  lmnD := lmnDivisors(lmn);
  QB := table(symmetric);
  for i from 1 to nops(lmnD) do
  for j from i to nops(lmnD) do
      QB[lmnD[i],lmnD[j]] := QuadBase(Pi/lmnD[i],Pi/lmnD[j]);
      od;
      od; 
  RETURN(QB);
  end:

SetMatrices := proc(lmn)
   local lamunu,ABC,Mats;
   lamunu := lammunu(lmn);
   ABC := FlmnTriple(lmn);
   Mats := table();
   Mats[_A] :=ABC[1];
   Mats[_iA] := inverse(Mats[_A]);
   Mats[_B] :=ABC[2];
   Mats[_iB] := inverse(Mats[_B]);
   Mats[_C] :=ABC[3];
   Mats[_iC] := inverse(Mats[_C]);
   Mats[_Ala] := evalm(ABC[1]^lamunu[1]);
   Mats[_iAla] := inverse(Mats[_Ala]);
   Mats[_Bmu] := evalm(ABC[2]^lamunu[2]);
   Mats[_iBmu] := inverse(Mats[_Bmu]);
   Mats[_Cnu] := evalm(ABC[3]^lamunu[3]);
   Mats[_iCnu] := inverse(Mats[_Cnu]);
   RETURN(Mats);
   end:

SetSchwartzEnv := proc(lmn,dig)
   global _glmn, _gTSides, _glamunu,_gMatrices,_gQBase,_glmndivs,_gDigits,_gGeoeps;
   local olddig;
   olddig := Digits;
   Digits := dig;
   _glmn := lmn;
   _glmndivs := lmnDivisors(lmn);
   _glamunu := lammunu(lmn);
   _gTSides := SetTSides(lmn); 
   _gMatrices := SetMatrices(lmn);
   _gQBase := SetQuadBases(lmn);
   _gDigits := dig; 
   _gGeoeps := 10.0^(5-dig);
    RETURN(NULL);
   end:

print(`done `);