with(linalg);

ortho:= proc(g,FN)
local o,i,ortho;
o[4]:=arccosh(((cosh(FN[3]/2)+cosh(FN[1]/2)*cosh(FN[2]/2))/sinh(FN[1]/2)/sinh(FN[2]/2)));
o[6]:=arccosh(((cosh(FN[2]/2)+cosh(FN[1]/2)*cosh(FN[3]/2))/sinh(FN[1]/2)/sinh(FN[3]/2)));
o[5]:=arccosh(((cosh(FN[1]/2)+cosh(FN[2]/2)*cosh(FN[3]/2))/sinh(FN[2]/2)/sinh(FN[3]/2)));
o[1]:=arccosh(sinh(o[4])^2*cosh(FN[2])-cosh(o[4])^2);
o[2]:=arccosh(sinh(o[4])^2*cosh(FN[1])-cosh(o[4])^2);
o[3]:=arccosh(sinh(o[6])^2*cosh(FN[1])-cosh(o[6])^2);
if g=2 then 
   for i from 7 to 12 do
       o[i]:=o[i-6];
   od;
fi;
if g>=3 then 
   #Q1
   o[8]:=arccosh(((cosh(FN[4]/2)+cosh(FN[1]/2)*cosh(FN[2]/2))/sinh(FN[1]/2)/sinh(FN[2]/2)));
   o[9]:=arccosh(((cosh(FN[2]/2)+cosh(FN[1]/2)*cosh(FN[4]/2))/sinh(FN[4]/2)/sinh(FN[1]/2)));
   o[11]:=arccosh(((cosh(FN[1]/2)+cosh(FN[2]/2)*cosh(FN[4]/2))/sinh(FN[2]/2)/sinh(FN[4]/2)));
   o[7]:=arccosh(sinh(o[8])^2*cosh(FN[2])-cosh(o[8])^2);
   o[10]:=arccosh(sinh(o[8])^2*cosh(FN[1])-cosh(o[8])^2);
   o[12]:=arccosh(sinh(o[9])^2*cosh(FN[1])-cosh(o[9])^2);
fi;
if g=3 then
   #P2
   o[14]:=arccosh(((cosh(FN[6]/2)+cosh(FN[3]/2)*cosh(FN[5]/2))/sinh(FN[3]/2)/sinh(FN[5]/2)));
   o[15]:=arccosh(((cosh(FN[5]/2)+cosh(FN[3]/2)*cosh(FN[6]/2))/sinh(FN[3]/2)/sinh(FN[6]/2)));
   o[17]:=arccosh(((cosh(FN[3]/2)+cosh(FN[5]/2)*cosh(FN[6]/2))/sinh(FN[5]/2)/sinh(FN[6]/2)));
   o[13]:=arccosh(sinh(o[14])^2*cosh(FN[5])-cosh(o[14])^2);
   o[16]:=arccosh(sinh(o[14])^2*cosh(FN[3])-cosh(o[14])^2);
   o[18]:=arccosh(sinh(o[15])^2*cosh(FN[3])-cosh(o[15])^2);

   #Q2
   o[20]:=arccosh(((cosh(FN[6]/2)+cosh(FN[4]/2)*cosh(FN[5]/2))/sinh(FN[4]/2)/sinh(FN[5]/2)));
   o[21]:=arccosh(((cosh(FN[5]/2)+cosh(FN[4]/2)*cosh(FN[6]/2))/sinh(FN[4]/2)/sinh(FN[6]/2)));
   o[23]:=arccosh(((cosh(FN[4]/2)+cosh(FN[5]/2)*cosh(FN[6]/2))/sinh(FN[5]/2)/sinh(FN[6]/2)));
   o[19]:=arccosh(sinh(o[20])^2*cosh(FN[5])-cosh(o[20])^2);
   o[22]:=arccosh(sinh(o[20])^2*cosh(FN[4])-cosh(o[20])^2);
   o[24]:=arccosh(sinh(o[21])^2*cosh(FN[4])-cosh(o[21])^2);
fi;
if g>3 then
   for i from 2 to g-2 do

       #Pi
       o[12*(i-1)+2]:=arccosh((cosh(FN[2*g-2+i]/2)+cosh(FN[i+1]/2)*cosh(FN[i+2]/2))/sinh(FN[i+1]/2)/sinh(FN[i+2]));
       o[12*(i-1)+3]:=arccosh((cosh(FN[2+i]/2)+cosh(FN[i+1]/2)*cosh(FN[i-2+2*g]/2))/sinh(FN[i+1]/2)/sinh(FN[i-2+2*g]));
       o[12*(i-1)+5]:=arccosh((cosh(FN[1+i]/2)+cosh(FN[i+2]/2)*cosh(FN[i-2+2*g]/2))/sinh(FN[i+2]/2)/sinh(FN[i-2+2*g]));
       o[12*(i-1)+1]:=arccosh((sinh(o[12*(i-1)+2])^2*cosh(FN[2+i])-cosh(o[12*(i-1)+2])^2));
       o[12*(i-1)+4]:=arccosh((sinh(o[12*(i-1)+2])^2*cosh(FN[1+i])-cosh(o[12*(i-1)+2])^2));
       o[12*(i-1)+6]:=arccosh((sinh(o[12*(i-1)+3])^2*cosh(FN[1+i])-cosh(o[12*(i-1)+3])^2));

       #Qi
              o[12*(i-1)+8]:=arccosh((cosh(FN[2*g-2+i]/2)+cosh(FN[i-1+g]/2)*cosh(FN[i+g]/2))/sinh(FN[i-1+g]/2)/sinh(FN[i+g]));
      o[12*(i-1)+9]:=arccosh((cosh(FN[i+g]/2)+cosh(FN[i-1+g]/2)*cosh(FN[i-2+2*g]/2))/sinh(FN[i-1+g]/2)/sinh(FN[i-2+2*g]));
       o[12*(i-1)+11]:=arccosh((cosh(FN[-1+i+g]/2)+cosh(FN[i+g]/2)*cosh(FN[i-2+2*g]/2))/sinh(FN[i+g]/2)/sinh(FN[i-2+2*g]));
       o[12*(i-1)+7]:=arccosh((sinh(o[12*(i-1)+8])^2*cosh(FN[i+g])-cosh(o[12*(i-1)+8])^2));
       o[12*(i-1)+10]:=arccosh((sinh(o[12*(i-1)+8])^2*cosh(FN[-1+i+g])-cosh(o[12*(i-1)+8])^2));
       o[12*(i-1)+12]:=arccosh((sinh(o[12*(i-1)+9])^2*cosh(FN[-1+i+g])-cosh(o[12*(i-1)+9])^2));
   od;
   
   #P(g-1)
   o[12*(g-2)+2]:=arccosh((cosh(FN[3*g-3]/2)+cosh(FN[3*g-4]/2)*cosh(FN[g]/2))/sinh(FN[3*g-4]/2)/sinh(FN[g+2]));
   o[12*(g-2)+3]:=arccosh((cosh(FN[3*g-4]/2)+cosh(FN[2+g]/2)*cosh(FN[3*g-3]/2))/sinh(FN[g]/2)/sinh(FN[3*g-3]));
   o[12*(g-2)+5]:=arccosh((cosh(FN[g]/2)+cosh(FN[3*g-3]/2)*cosh(FN[3*g-4]/2))/sinh(FN[3*g-4]/2)/sinh(FN[3*g-3]));
   o[12*(g-2)+1]:=arccosh((sinh(o[12*(g-2)+2])^2*cosh(FN[3*g-4])-cosh(o[12*(g-2)+2])^2));
   o[12*(g-2)+4]:=arccosh((sinh(o[12*(g-2)+2])^2*cosh(FN[g])-cosh(o[12*(g-2)+2])^2));
   o[12*(g-2)+6]:=arccosh((sinh(o[12*(g-2)+3])^2*cosh(FN[g])-cosh(o[12*(g-2)+3])^2));

 #Q(g-1)
   o[12*(g-2)+8]:=arccosh((cosh(FN[3*g-3]/2)+cosh(FN[3*g-4]/2)*cosh(FN[2*g-2]/2))/sinh(FN[3*g-4]/2)/sinh(FN[2*g-2]));
   o[12*(g-2)+9]:=arccosh((cosh(FN[3*g-4]/2)+cosh(FN[-2+2*g]/2)*cosh(FN[3*g-3]/2))/sinh(FN[-2+2*g]/2)/sinh(FN[3*g-3]));
   o[12*(g-2)+11]:=arccosh((cosh(FN[-2+2*g]/2)+cosh(FN[3*g-3]/2)*cosh(FN[3*g-4]/2))/sinh(FN[3*g-4]/2)/sinh(FN[3*g-3]));
   o[12*(g-2)+7]:=arccosh((sinh(o[12*(g-2)+8])^2*cosh(FN[3*g-4])-cosh(o[12*(g-2)+8])^2));
   o[12*(g-2)+10]:=arccosh((sinh(o[12*(g-2)+8])^2*cosh(FN[-2+2*g])-cosh(o[12*(g-2)+8])^2));
   o[12*(g-2)+12]:=arccosh((sinh(o[12*(g-2)+9])^2*cosh(FN[-2+2*g])-cosh(o[12*(g-2)+9])^2));
fi;
RETURN(o);
end;



Nnum:=proc(g,FN)
local o,y,i,t,N;
o:=convert(ortho(g,FN),array);
y:=FN[1];
for i from 2 to 3*g-3 do
    if evalf(y)>evalf(FN[i]) then y:=FN[i] fi;
od;
t:=o[1];
for i from 2 to 12*g-12 do
    if evalf(t)>evalf(o[i]) then t:=o[i] fi; 
od;
RETURN(y/t);
end;




N:=proc(g,FN)
local o,y,i,t,N;
o:=convert(ortho(g,FN),array);
y:=FN[1];
for i from 2 to 3*g-3 do
    if signum(y-FN[i])=1 then y:=FN[i] fi;
od;
t:=o[1];
for i from 2 to 12*g-12 do
    if signum(t-o[i])=1 then t:=o[i] fi; 
od;
RETURN(simplify(round(evalf(simplify(y/t)))));
end;

#Recherche des chemins de longueur n
paths:=proc(genre,n)
local Q,i,P,c1,c,k;
if evalb(n=1) then 
if genre=2 then
Q:={[[1,'h']],[[1,'b']],[[2,'h']],[[2,'b']],[[3,'h']],[[3,'b']]};
fi;
if genre>2 then      Q:={[[1,'h']],[[1,'b']],[[2,'h']],[[2,'b']],[[3*genre-3,'h']],[[3*genre-3,'b']],[[3*genre-4,'h']],[[3*genre-4,'b']]};
fi;
for i from 3 to genre do
    Q:=Q union {[[i,'d']],[[i,'g']]};
od;
for i from genre+1 to 2*genre-2 do
    Q:=Q union {[[i,'d']],[[i,'g']]};
od;
for i from 2*genre-1 to 3*genre-5 do
    Q:=Q union {[[i,'h']],[[i,'b']]};
od;
fi;
if evalf(signum(n-1))=1. then
   P:=paths(genre,n-1);
   Q:=P;
   c1:=array(1..n);
   for c in P do
       if nops(c)=n-1 then
       for k from 1 to n-1 do
           c1[k]:=c[k];
       od; 
       if genre=2 then
          if (c[n-1]=[1,'h']) or c[n-1]=[2,'h'] or c[n-1]=[3,'h'] then
             c1[n]:=[1,'b'];
             Q:=Q union {convert(c1,list)};
             c1[n]:=[2,'b'];
             Q:=Q union {convert(c1,list)};
             c1[n]:=[3,'b'];
             Q:=Q union {convert(c1,list)};
          fi;
          if c[n-1]=[1,'b'] or c[n-1]=[2,'b'] or c[n-1]=[3,'b']then
             c1[n]:=[1,'h'];
             Q:=Q union {convert(c1,list)};
             c1[n]:=[2,'h'];
             Q:=Q union {convert(c1,list)};
             c1[n]:=[3,'h'];
             Q:=Q union {convert(c1,list)};
          fi
        fi;
        if genre=3 then
           if (c[n-1]=[1,'h']) or c[n-1]=[2,'h'] or c[n-1]=[3,'g'] then
             c1[n]:=[1,'b'];
             Q:=Q union {convert(c1,list)};
             c1[n]:=[2,'b'];
             Q:=Q union {convert(c1,list)};
             c1[n]:=[3,'d'];
             Q:=Q union {convert(c1,list)};
           fi; 
           if (c[n-1]=[5,'h']) or c[n-1]=[6,'h'] or c[n-1]=[3,'d'] then
             c1[n]:=[5,'b'];
             Q:=Q union {convert(c1,list)};
             c1[n]:=[6,'b'];
             Q:=Q union {convert(c1,list)};
             c1[n]:=[3,'g'];
             Q:=Q union {convert(c1,list)};
           fi;
           if (c[n-1]=[1,'b']) or c[n-1]=[2,'b'] or c[n-1]=[4,'g'] then
             c1[n]:=[1,'h'];
             Q:=Q union {convert(c1,list)};
             c1[n]:=[2,'h'];
             Q:=Q union {convert(c1,list)};
             c1[n]:=[4,'d'];
             Q:=Q union {convert(c1,list)};
           fi;
           if (c[n-1]=[5,'b']) or c[n-1]=[6,'b'] or c[n-1]=[4,'d'] then
             c1[n]:=[5,'h'];
             Q:=Q union {convert(c1,list)};
             c1[n]:=[6,'h'];
             Q:=Q union {convert(c1,list)};
             c1[n]:=[4,'g'];
             Q:=Q union {convert(c1,list)};
           fi;
        fi;
        if genre>3 then
           #P1 
           if (c[n-1]=[1,'h']) or c[n-1]=[2,'h'] or c[n-1]=[3,'g'] then
             c1[n]:=[1,'b'];
             Q:=Q union {convert(c1,list)};
             c1[n]:=[2,'b'];
             Q:=Q union {convert(c1,list)};
             c1[n]:=[3,'d'];
             Q:=Q union {convert(c1,list)};
           fi; 
           #Q1 
           if (c[n-1]=[1,'b']) or c[n-1]=[2,'b'] or c[n-1]=[genre+1,'g'] then
             c1[n]:=[1,'h'];
             Q:=Q union {convert(c1,list)};
             c1[n]:=[2,'h'];
             Q:=Q union {convert(c1,list)};
             c1[n]:=[genre+1,'d'];
             Q:=Q union {convert(c1,list)};
           fi ;
           #Pi
           for i from 2 to genre-2 do
             if (c[n-1]=[i+1,'d']) or c[n-1]=[i+2,'g'] or c[n-1]=[2*genre-3+i,'h'] then
               c1[n]:=[i+1,'g'];
               Q:=Q union {convert(c1,list)};
               c1[n]:=[i+2,'d'];
               Q:=Q union {convert(c1,list)};
               c1[n]:=[2*genre-3+i,'b'];
               Q:=Q union {convert(c1,list)};
             fi; 
           od;
           #Qi
           for i from 2 to genre-2 do
             if (c[n-1]=[genre+i-1,'d']) or c[n-1]=[genre+i,'g'] or c[n-1]=[2*genre-3+i,'b'] then
               c1[n]:=[genre+i-1,'g'];
               Q:=Q union {convert(c1,list)};
               c1[n]:=[genre+i,'d'];
               Q:=Q union {convert(c1,list)};
               c1[n]:=[2*genre-3+i,'h'];
               Q:=Q union {convert(c1,list)};
             fi; 
           od;
           #P(g-1)
         
             if (c[n-1]=[3*genre-3,'h']) or c[n-1]=[3*genre-4,'h'] or c[n-1]=[genre,'d'] then
               c1[n]:=[3*genre-3,'b'];
               Q:=Q union {convert(c1,list)};
               c1[n]:=[3*genre-4,'b'];
               Q:=Q union {convert(c1,list)};
               c1[n]:=[genre,'g'];
               Q:=Q union {convert(c1,list)};
             fi ;
           #Q(g-1)
           if (c[n-1]=[3*genre-3,'b']) or c[n-1]=[3*genre-4,'b'] or c[n-1]=[2*genre-2,'d'] then
               c1[n]:=[3*genre-3,'h'];
               Q:=Q union {convert(c1,list)};
               c1[n]:=[3*genre-4,'h'];
               Q:=Q union {convert(c1,list)};
               c1[n]:=[2*genre-2,'g'];
               Q:=Q union {convert(c1,list)};
             fi ;
     fi;
  fi;
  od;
fi;
RETURN(Q);
end;

 
loops:=proc(genre,n)
local Q,p;
Q:=convert(paths(genre,n),set);
for p in Q do
    if nops(p)=1 then Q:=Q minus {p};fi;
    if p[1]<>p[nops(p)] then Q:=Q minus {p};fi;od;
RETURN(Q);
end;

get_poss:=proc(a,b,c,genre,F)
local Q;
if genre=2 then
   if ((b[1]<>a[1]) and (a[1]=c[1])) then 
      if F[3+b[1]]<>0 then 
         Q:={[signum(F[3+b[1]])*Pi/2,abs(F[3+b[1]])],[-signum(F[3+b[1]])*Pi/2,F[b[1]]-abs(F[3+b[1]])]};
      else
         Q:={[Pi/2,0],[Pi/2,F[b[1]]],[-Pi/2,F[b[1]]]};
      fi;
   fi;
   if (b[1]<>a[1] and a[1]<>c[1] and b[1]<>c[1]) then 
      if F[3+b[1]]<>0 then 
         Q:={[signum(F[3+b[1]])*Pi/2,F[b[1]]/2+abs(F[3+b[1]])],[-signum(F[3+b[1]])*Pi/2,F[b[1]]/2-abs(F[3+b[1]])]};
      else 
         Q:={[Pi/2,F[b[1]]/2],[-Pi/2,F[b[1]]/2]};
      fi;
   fi;

   if ((b[1]=a[1] or b[1]=c[1]) and (a[1]<>c[1])) then 
      if signum(-F[b[1]]/2-F[3+b[1]])=-1 and signum(F[3+b[1]]+F[b[1]]/4)=-1 then 
         Q:={[Pi/2,3*F[b[1]]/4+F[3+b[1]]],[Pi/2,5*F[b[1]]/4+F[3+b[1]]],[-Pi/2,-F[b[1]]/4-F[3+b[1]]],[-Pi/2,F[b[1]]/4-F[3+b[1]]]};
      fi;
      if signum(-F[b[1]]/4-F[3+b[1]])=-1 and signum(F[3+b[1]]-F[b[1]]/4)=-1 then
         Q:={[Pi/2,F[b[1]]/4+F[3+b[1]]],[Pi/2,3*F[b[1]]/4+F[3+b[1]]],[-Pi/2,3*F[b[1]]/4-F[3+b[1]]],[-Pi/2,F[b[1]]/4-F[3+b[1]]]};
      fi;
      if signum(F[b[1]]/4-F[3+b[1]])=-1 and signum(F[3+b[1]]-F[b[1]]/2)=-1 then
         Q:={[Pi/2,F[b[1]]/4+F[3+b[1]]],[Pi/2,-F[b[1]]/4+F[3+b[1]]],[-Pi/2,3*F[b[1]]/4-F[3+b[1]]],[-Pi/2,5*F[b[1]]/4-F[3+b[1]]]};
      fi;
      if F[3+b[1]]= -F[b[1]]/4 then
         Q:={[Pi/2,F[b[1]]/2],[Pi/2,F[b[1]]],[-Pi/2,0],[-Pi/2,F[b[1]]/2],[Pi/2,0],[-Pi/2,F[b[1]]]};
      fi;
      if F[3+b[1]]= F[b[1]]/4 then
         Q:={[Pi/2,F[b[1]]/2],[Pi/2,F[b[1]]],[-Pi/2,0],[-Pi/2,F[b[1]]/2],[Pi/2,0],[-Pi/2,F[b[1]]]};
      fi;
      if F[3+b[1]]= F[b[1]]/2 then
         Q:={[Pi/2,3*F[b[1]]/4],[Pi/2,F[b[1]]/4],[-Pi/2,F[b[1]]/4],[-Pi/2,3*F[b[1]]/4]};
      fi;
   fi;

   if (a[1]=b[1] and b[1]=c[1]) then 
      if  F[3+b[1]]<>0 then 
Q:={[signum(F[3+b[1]])*Pi/2,abs(F[3+b[1]])],[signum(F[3+b[1]])*Pi/2,F[1]/2+abs(F[3+b[1]])],[-signum(F[3+b[1]])*Pi/2,F[1]/2-abs(F[3+b[1]])],[-signum(F[3+b[1]])*Pi/2,F[1]-abs(F[3+b[1]])]};
      else
Q:={[Pi/2,0],[Pi/2,F[1]/2],[-Pi/2,F[1]/2],[-Pi/2,F[1]],[Pi/2,F[1]]};
      fi;
   fi;
fi;
if genre>2 then
   print("Non implemente");
   Q:={};
fi;
RETURN(Q);
end;



get_val:=proc(i,L,N,genre,F,T,Q) local R,P;
R:=Q;
if i=1 then 
   for P in get_poss(L[N-1],L[1],L[2],genre,F) do 
       T[1]:=P;
       R:=get_val(2,L,N,genre,F,T,R);
   od;
fi;
if 1<i and i<=N-1 then 
   for P in get_poss(L[i-1],L[i],L[i+1],genre,F) do 
       T[i]:=P;
       R:=get_val(i+1,L,N,genre,F,T,R);
   od;
fi;
if i=N then
   R:=R union {convert(T,array)};   
fi;
RETURN(R);
end;



get_class:=proc(F)
local genre,i,R,L,l,T;
genre:=1/6*(nops(F)+6);
i:=eval(N(genre,F));
R:={};
L:=loops(genre,i+1);
for l in L do
   
   T:=array(1..nops(l)-1);
   R:=R union {[get_val(1,l,nops(l),genre,F,T,{}),l]};
od;
RETURN(R);

end;

R:=theta->matrix(2,2,[exp(I*theta/2),0,0,exp(-I*theta/2)]);
T:=t->matrix(2,2,[cosh(t/2),sinh(t/2),sinh(t/2),cosh(t/2)]);

get_syst:=proc(F)
local epsilon,genre,c,o,s,S,i,C,val,M;
#DEBUG();
epsilon:=10^(-6);
genre:=1/6*(nops(F)+6);

C:=get_class(F);
o:=ortho(genre,F);
s:=F[1];
S:={};

for i from 1 to 3*genre-3 do
    if evalf(F[i]- s+epsilon)<0 then
       S:={[i,'h']};
       s:=F[i];
    fi;
    if evalf(abs(F[i]- s)-epsilon)<=0 then
       S:=S union {[i,'h']};
    fi;
od;

for c in C do
     #DEBUG(S,evalf(s),c[2]); 
     for val in c[1] do          
        #DEBUG("val in c[1]");
        M:=get_matrix(val,c[2],o);
        #DEBUG(c[2]);
        t:=(2*arccosh(trace(evalm(M))/2));
        tnum:=evalf(t);
        #DEBUG(evalf(t),evalf(s),c[2]);

        if Re((evalf(abs(tnum-evalf(s))-epsilon)))<=0 then
           #DEBUG("one another");
           S:=S union {[convert(val,array),c[2]]};
        fi;
        if Re((evalf((abs(tnum)-abs(evalf(s))+epsilon))))<0 then 
           #DEBUG("new one");     
           S:={[convert(val,array),c[2]]};
           s:=t;
        fi;
    od;
od;
return([s,S]);
end;

get_matrix:=proc(val,loop,o)
local M,i;
    M:=matrix(2,2,[1,0,0,1]);
    for i from 1 to nops(loop)-1 do
        if loop[i+1][1]=loop[i][1] then
           M:=evalm(R(-val[i][1])&*T(val[i][2])&*R(val[i][1])&*T(o[i])&*M);   
        fi;
        if (loop[i][1]=1 and loop[i+1][1]=2) or (loop[i][1]=2 and loop[i+1][1]=1) then
           M:=evalm(R(-val[i][1])&*T(val[i][2])&*R(val[i][1])&*T(o[4])&*M);
        fi;
        if (loop[i][1]=1 and loop[i+1][1]=3) or (loop[i][1]=3 and loop[i+1][1]=1) then
           M:=evalm(R(-val[i][1])&*T(val[i][2])&*R(val[i][1])&*T(o[6])&*M);
        fi;
        if (loop[i][1]=2 and loop[i+1][1]=3) or (loop[i][1]=3 and loop[i+1][1]=2) then
           M:=evalm(R(-val[i][1])&*T(val[i][2])&*R(val[i][1])&*T(o[5])&*M);
        fi;
   od;
   RETURN(M);

end;