with(plots):
Digits:=400;

## fonction de calcul de la suite de Sturm-Habicht

StHa:= proc(a,b,x)
  local suc,p,q,n1,n2,n,m,lp,lq,r,R,S0,S1,d,NS,Sr,j;
  p:=expand(a); q:=expand(b);
  if q=0 then RETURN(x,p); fi;
  n1:=degree(p,x); n2:=degree(q,x); lp:=lcoeff(p,x); lq:=lcoeff(q,x); 
  suc:=p,q;
  if n1>n2 then
    n:=n1-1; d:=n-n2;
      if d<>0 then 
        S1:=((-1)**(d*(d+1)/2))*(lp*lq)**d*q; S0:=((-1)**((d+2)*(d+3)/2))*lp**d*prem(p,q,x);
        if S0<>0 then suc:=suc,S1,S0 else suc:=suc,S1 fi;
      else S1:=q; S0:= -prem(p,q,x); if S0<>0 then suc:= suc,S0 fi;
      fi;
   j:=n2-1;
else 
   n:=n2; S1:=q; S0:= -rem(lq**2*p,q,x); j:=n-1; if S0<>0 then suc:=suc,S0 fi;
fi;
R:=lcoeff(S1,x); if S0=0 then r:=-1 else r:=degree(S0,x) fi;
while r>=0 do
   if r>0 then
      if r=j then 
         NS:='NS'; divide(prem(S1,S0,x),((-R)**(j-r+2)),NS); NS:=expand(NS)*((-1)**((j-r+2)*(j-r+3)/2));
         if NS<>0 then suc:=suc,NS fi;
         S1:=S0; S0:=NS; R:=lcoeff(S1,x); j:=r-1;
      else 
         for m from r+1 to j-1 do suc:=suc,0 od: Sr:='Sr';
         divide(prem(S1,S0,x),(R**(j-r+2)),Sr); S1:='S1';
         divide(((lcoeff(S0,x)**(j-r))*S0),(R**(j-r)),S1);
         S1:=S1*((-1)**((j-r)*(j-r+1)/2)); S0:=Sr*((-1)**((j-r+2)*(j-r+3)/2));
         if S0<>0 then suc:=suc,S1,S0 else suc:=suc,S1 fi;
         R:=lcoeff(S1,x); j:= r-1;
     fi
   else 
      if j>0 then
         S1:='S1'; divide(((lcoeff(S0,x)**j)*S0),(R**j),S1);
         S1:=S1*((-1)**(j*(j+1)/2)); suc:=suc,S1;
      fi:
      S0:=0;
   fi:
   if S0<>0 then r:=degree(S0,x) else r:=-1 fi:
od:
suc;
end:

## calcul des listes des coefficients sous-résultants et des gamma
## F un polynôme en les variables x et y
## Cette fonction retourne sr qui est la liste des listes des coefficients 
## sous-résultants et gam est la liste des polynômes gamma.
Gam:=proc(F)
  local Phi,gam,sr,L,n,i,j,k,m;
  L:=[StHa(F,diff(F,y),y)];
  n:= nops(L);
  for i from 1 by 1 to n do 
    m[i]:=degree(L[i],y);
    for j from 0 by 1 to m[i] do
      sr[n-i][j]:=coeff(L[i],y,j);
    end do;
  end do;
  Phi[0]:=quo(sr[0][0],gcd(sr[0][0],diff(sr[0][0],x)),x);
  for j from 1 by 1 to n-1 do
    Phi[j]:=gcd(Phi[j-1],sr[j][j]);
    gam[j]:=quo(Phi[j-1],Phi[j],x);
  end do;
gam:=convert(gam,'list');
  RETURN([sr,gam]);
end proc:


# generic test
IsGeneric := proc(sr,gam)
  local k, r, n, cond, i;
  n := nops(gam);
  k := 1;
  cond := true;
  while ((k < n) and (cond = true)) do
    i := 0;
    if degree(gam[k]) > 0 then
      while ((i < k) and (cond = true)) do
	 r[k][i]:=rem(k*(k-i)*(sr[k][i])*(sr[k][k])-(i+1)*(sr[k][k-1])*(sr[k][i+1]),gam[k],x);
print(r[k][i]);
         if ( r[k][i] <> 0 ) then cond:=false else i:=i+1
	 fi;
      od;
    fi;
k:=k+1;
  od;
  RETURN(cond);
end:

## Calcul des paramétrisations rationnelles associées à chaque gamma
## sr est liste des listes des coefficients sous-résultants et gam 
## est la liste des polynômes gamma
## Cette fonction retourne beta qui est la liste des parametrisations rationnelles
## associées au gamma. Aisni beta[i] est associée à gam[i].
ComputeBetaCrit := proc(sr,gam)
  local i, beta;
  for i from 1 to nops(gam)do
 if gam[i]<>1 then
    beta[i] := -(sr[i][i-1])/(i*sr[i][i]);
	 else beta[i]:=nini;
 fi;
od;
  RETURN(beta);
end:

## Calcul des point critique
## gam est la liste des polynômes gamma de la factorisation du résultant
## beta et la paramétrisation rationnelle associée à chaque gamma
## ATTENTION ses listes sont ordonnées
## retourne une liste de deux listes : la première liste est la liste des valeurs 
## critiques et la deuxième est la liste des ordonnées correspondantes dans le 
## même ordre
##INUTILE
CriticalPoints := proc(gam,beta)
  local pc,i,co;
  pc := [seq([fsolve(gam[i])],i=1..nops(gam))];
  co := [seq([seq(eval(beta[i],x=pc[i][j]),j=1..nops(pc[i]))],i=1..nops(gam))];
  RETURN([[seq(op(pc[i]),i=1..nops(pc))],[seq(op(co[i]),i=1..nops(co))]]);
end:

InsertInFib := proc(a,L)
  local i,j;
  if L = [] then
    RETURN([a]);
  else
    i := 1;
	while ( (L[i] > a) and (i < nops(L)+1) ) do
	  i := i+1;
	od;
	if ( i = nops(L) ) then
      RETURN([op(L),a]);
	else
	  RETURN([seq(L[j],j=1..i-1),a,seq(L[j],j=i..nops(L))]);
	fi;
  fi;
end:

deflate := proc(F,gamma,beta,i)
  local f,v,d;
d:=degree(F,y);
  f := rem(F,gamma,x);
v:=expand(subs(y=(y+beta),F));
v:=numer(v);
v:=quo(v,y^(i+1),y);
  f :=rem(v,gamma,x);
  RETURN(f);
end:

CritFib := proc(f,a,b)
  local L, Up, Down, j,i;
  L := [fsolve(subs(x=a,f),y)];
  Up := [];
  Down := [];
  for j from 1 to nops(L) do
    if L[j] > 0 then 
	  Up := InsertInFib(L[j],Up);
	else
	  Down := InsertInFib(L[j],Down);
	fi;
  od;
for i from 1 to nops(Up) do Up[i]:=Up[i]+b;
od;
for i from 1 to nops(Down) do Down[i]:=Down[i]+b;
od;
  RETURN([Up, b, Down]);
end:

CriticalFibers := proc(F,gam,beta)
  local pc, i, j, fib, fibs, f, co;
  pc := [seq([fsolve(gam[i])],i=1..nops(gam))];
  co := [seq([seq(eval(beta[i],x=pc[i][j]),j=1..nops(pc[i]))],i=1..nops(gam))];
  fibs := [];
  for i from 1 to nops(pc) do
    f := deflate(F,gam[i],beta[i],i);
	fib := [];
	for j from 1 to nops(pc[i]) do
	  fib := [op(fib),CritFib(f,pc[i][j],co[i][j])];
	od;
	fibs := [op(fibs),fib];
  od;
  RETURN([[seq(op(pc[i]),i=1..nops(pc))],[seq(op(fibs[i]),i=1..nops(fibs))]]);
end:

## recherche du plus grand élément dans une liste.
mon_sup := proc(L,i,j)
  local r,k,ind;
  r := L[i];
  ind := i;
  for k from i+1 to j do
    if (L[k] > r) then
      r := L[k];
      ind := k;
    fi:
  od;
  RETURN(ind);
end:

## Calcul des ordonnées des points d'une fibre régulière
ComputeRegFiber := proc(F,v)
RETURN(sort([fsolve(subs(x=v,F),y)], `>`));
end:

#trier une liste par ordre décroissant
## L est une liste de deux listes. La fonction trie la première liste par permuttation
## et reporte les opérations sur la deuxième liste
MonTri := proc(L)
  local A,B,a,i,j,v;
  A := L[1];
  B := L[2];
  a := nops(A);
  for i from 1 to a-1 do
    j := mon_sup(A,i,a);  
    v := A[j];
    A[j] := A[i];
    A[i] := v;
    v := B[j];
    B[j] := B[i];
    B[i] := v;
  od;
  RETURN([A,B]);
end:

## connection d'une fibre critique à une fibre régulière
ConnectCritToReg := proc(C,R)
  local Up, reg, Down, d, r, i, j, u;
  Up:=sort(C[2], `>`);
  u := nops(Up);
print("u",u);
  Down:=sort(C[4], `>`);
  d := nops(Down);
  reg := R[2];
  r := nops(reg);
  j := r - d - u;
  RETURN([seq([ [C[1],Up[i]],[R[1],reg[i]] ],i=1..u),seq([ [C[1],C[3]],[R[1],reg[u+i]] ],i=1..j), seq([ [C[1],Down[d-i+1]],[R[1],reg[r-i+1]] ],i=1..d)]);
end:

RegValues := proc(cv)
  local res, i;
if nops(cv)>1 then
  res := [seq((cv[i+1]+cv[i])/2,i=1..(nops(cv)-1))];
  res := [cv[1]-(res[1]-cv[1]),op(res),cv[nops(cv)]+(cv[nops(cv)]-res[nops(cv)-1])];
else res:=[cv[1]-1,cv[1]+1];
fi;
  RETURN(res);
end:

MaTopo := proc(F)
  local G,gamma, sr, cf, beta, cv, complexe, rv, rf,C,R;
  #calcul des fibre critiques
  G := Gam(F);
  gamma := G[2];
  sr := G[1];
  beta := ComputeBetaCrit(sr,gamma);
  cf := CriticalFibers(F,gamma,beta);
  cf := MonTri(cf);
  cv := cf[1];
  cf := cf[2];
  ## calcul des fibres régulières
  rv := RegValues(cv);
  rf := [seq(ComputeRegFiber(F,rv[i]),i=1..nops(rv))];
  ## mise à jour des structures de données des fibres pour la connection
  C := [seq([cv[i],op(cf[i])], i = 1..nops(cv))];
  R := [seq([rv[i],rf[i]],i=1..nops(rf))];
  complexe := [seq(ConnectCritToReg(C[i],R[i+1]),i=1..nops(C)),seq(ConnectCritToReg(C[nops(C)-i+1],R[nops(C)-i+1]), i=1..nops(C))];
  RETURN(complexe); 
end:

Draw := proc(F)
  local comp,C,i;
  comp := MaTopo(F);
  C := [];
  for i from 1 to nops(comp) do
    if ( nops(comp[i]) > 0 ) then
      C  := [op(comp[i]), op(C)];
    fi;
  od:
  display(CURVES(op(C)));
end: