with(plots):
with(LinearAlgebra):
with(VectorCalculus):
Digits:=300:

## 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:



## La fonction SquareFree calcul la partie sans facteurs carrés d'un polynôme bivarié
## F est un polynôme en les variables x et y.

SquareFree:=proc(F)
  local m,i,f,P,W,R,V;

	m:=degree(F,y);

for i from 0 by 1 to m do 
  f[i]:=coeff(F,y,i);
od;
	P:=f[0];
	i:=1;
while (i <= m) do
	P:=gcd(P,f[i]);
  	i:=i+1;
od;

	divide(F,P,W);
	R:=gcd(W,diff(W,y));
	divide(W,R,V);
	P:=quo(P,gcd(P,diff(P,x)),x);

RETURN(P*V);
end:



## calcul des listes des coefficients sous-résultants et des deltas
## P et Q sont des polynômes en les variables x, y et z.
## Cette fonction retourne SR qui est la liste des listes des coefficients 
## sous-résultants et del qui est la liste des polynômes deltas.

Delt:=proc(P,Q)
  local Phi,del,SR,L,n,i,j,k;
	L:=[StHa(P,Q,z)];
	n:= nops(L);
  for i from 1 by 1 to n do 
    for j from 0 by 1 to (n-i) do
	SR[n-i][j]:=coeff(L[i],z,j);
    od;
  od;
	Phi[0]:=SquareFree(SR[0][0]);
  for j from 1 by 1 to n-1 do
	Phi[j]:=gcd(Phi[j-1],SR[j][j]);
	divide(Phi[j-1],Phi[j],del[j]);
  od;
	del:=convert(del,'list');
  RETURN([SR,del,Phi[0],L]);
end:


## La fonction IsPseudoGeneric test si une courbe 3D donnée est en position
## pseudo-generique. Elle prends en entrée les listes des coefficients
## sous-resultants et des polynômes deltas associées à notre courbe 3D. 


IsPseudoGeneric := proc(SR,del)
  local k, r, n, cond, i;
	n := nops(del);
 	k :=1;
 	cond := true;
  while ((k <= n) and (cond = true)) do
 	i := 0;
     if degree(del[k]) > 0 then
      while ((i < k) and (cond = true)) do
	cond:= divide((k*(k-i)*(SR[k][i])*(SR[k][k])-(i+1)*(SR[k][k-1])*(SR[k][i+1])),del[k]);
	i:=i+1;
      od;
     fi;
	k:=k+1;
  od;
  RETURN(cond);
end:




## calcul des listes des coefficients sous-résultants et des gammas
## 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;
	L:=[StHa(F,diff(F,y),y)];
	n:= nops(L);

  for i from 1 by 1 to n do 
    for j from 0 by 1 to (n-i) do
	sr[n-i][j]:=coeff(L[i],y,j);
    od;
  od;

	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);
  od;

	gam:=convert(gam,'list');
  RETURN([sr,gam]);
end:


## La fonction IsGeneric test si la courbe plane, projection de nôtre
## courbe 3D, est en position générique.

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],x) > 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);
       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 degree(gam[i],x)>0 then
	beta[i] := -rem((sr[i][i-1]),gam[i],x)/rem((i*sr[i][i]),gam[i],x);

   else
	beta[i]:=nini;
   fi;
  od;
	beta:=convert(beta,'list');
  RETURN(beta);
end:



## Distinction entre singularités et singularités apparentes
## et test de généricité de la position d'une courbe gauche réelle. 

IsGGeneric:= proc(SR,gam,beta,L)
  local gamm,h,SRR,R,RR,k,i,j,v,m,n,f,AppSing,g,w,F,FF,cond,lamda,betaa,LL,Appbeta,Applamda,b;
	
	k:=nops(beta);
	n:=nops(L);
for h from 1 by 1 to k do
   if degree(gam[h],x)>0 then 
	SRR[1][1]:=numer(subs(y=beta[h],SR[1][1]));
	v:=gcd(SRR[1][1],gam[h]);
	gamm[h][1]:=quo(gam[h],v,x);

   if degree(gamm[h][1],x)>0 then
	betaa[h][1]:=rem(numer(beta[h]),gamm[h][1],x)/rem((denom(beta[h])),gamm[h][1],x);
   else
	betaa[h][1]:=beta[h];
   fi;

   if degree(gamm[h][1],x)>0 then 
	lamda[h][1]:=numer(subs(y=betaa[h][1],L[n-1]));
	lamda[h][1]:=rem(lamda[h][1],gamm[h][1],x);
	lamda[h][1]:=-coeff(lamda[h][1],z,0)/coeff(lamda[h][1],z,1);
   fi;


	g[h]:=v;
  if degree(g[h],x)>0 then
	b[h]:=rem(numer(beta[h]),g[h],x)/rem(denom(beta[h]),g[h],x);
  else 
	b[h]:=beta[h];
  fi;


   else
	g[h]:=1;

	b[h]:=1;
   fi;
od;
for i from 1 by 1 to k do
	j:=2; 
    while(degree(g[i],x)>0) do
	SRR[j][j]:=numer(subs(y=b[i],SR[j][j]));
	f:=gcd(SRR[j][j],g[i]);
	v:=quo(g[i],f,x);
	w:=v;
	
  for m from 0 by 1 to (j-1) do
      	R[j][m]:=((j)*(j-m)*(SR[j][m])*(SR[j][j])-(m+1)*(SR[j][j-1])*(SR[j][m+1]));
      	RR[j][m]:=numer(subs(y=b[i],R[j][m]));
      	v:=gcd(RR[j][m],v);
  od;

   if (degree(v,x)>0) then
  	gamm[i][j]:=v;
	betaa[i][j]:=rem(numer(b[i]),gamm[i][j],x)/rem(denom(b[i]),gamm[i][j],x);
	LL[n-j]:=numer(subs(y=betaa[i][j],L[n-j]));
	lamda[i][j]:=-rem(coeff(LL[n-j],z,(j-1)),gamm[i][j],x)/rem(j*coeff(LL[n-j],z,j),gamm[i][j],x);
	w:=quo(w,v,x);
    else
	gamm[i][j]:=1;
	betaa[i][j]:=1;
	lamda[i][j]:=1;
	
   fi;


   if (degree(w,x)>0) then 
	AppSing[i][j]:=w;
	Appbeta[i][j]:=rem(numer(b[i]),AppSing[i][j],x)/rem(denom(b[i]),AppSing[i][j],x);
	Applamda[i][j]:=numer(subs(y=Appbeta[i][j],L[n-j]));
	Applamda[i][j]:=rem(Applamda[i][j],AppSing[i][j],x);
   else
	AppSing[i][j]:=1;
	Appbeta[i][j]:=1;
	Applamda[i][j]:=1; 
   fi;
  	j:=j+1;
  	g[i]:=f;
   od; 
	i:=i+1;
od;

	cond := true;
    for i from 1 by 1 to k do
	j:=2; 
    while((j < n) and (cond = true)) do
      if (degree(AppSing[i][j],x) > 0) then
	F:=resultant(L[n-j],diff(L[n-j],z),z);
        FF:=numer(subs(y=Appbeta[i][j],F));
        FF:=gcd(AppSing[i][j],FF);

       if ( degree(FF,x)> 0 ) then
	cond:=false;
       else
	j:=j+1;
       fi;
      else
	j:=j+1; 
      fi; 
  od;
od;

RETURN([cond,AppSing,Appbeta,Applamda,gamm,betaa,lamda]);    
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:



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:


## 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:


##trier une liste par ordre décroissant
## L est une liste de deux listes. La fonction trie la première liste par permutation
## 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:


CritFib := proc(f,a,b,c)
  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,c], Down]);
end:

## cette fonction permet de relever les points d'une fibre x-critique plane dans l'espace.

LiftCritFib:=proc(fib,a,SR)
 local Up,Down,b,i;

	Up:=fib[1];
	Down:=fib[3];
	b:=fib[2];
	Up:=[[op(Up)],[seq(eval(eval(-SR[1][0]/SR[1][1],x=a),y=Up[i]),i=1..nops(Up))]];
	Down:=[[op(Down)],[seq(eval(eval(-SR[1][0]/SR[1][1],x=a),y=Down[i]),i=1..nops(Down))]];
	Up:=MonTri(Up);
	Down:=MonTri(Down);

 RETURN([Up,b,Down]);
end:

## Cette fonction permet de trier les points d'une fibre x-critique en apparence.
## fib est la fibre x-critique en apparence
## a est la valeur x-critique
## Jac est le noyeau de la matrice jacobienne associé à P et Q.
AppSingTri:=proc(a,Jac,fib)
 local f,g,b,i,L;
	g:=fib[2];
	b:=fib[1];
	L:=[];

  for i from 1 by 1 to nops(g) do
	f[i]:=eval(eval(eval(Jac,x=a),y=b),z=g[i]);
  if (f[i][1] < 0) then
	f[i]:=-f[i];
  fi;
	f[i]:=f[i]/sqrt((f[i]).(f[i]));
	L:=[op(L),f[i][2]];
  od;
	L:=MonTri([L,g]);
RETURN([b,L[2]]);
end:






## Calcul des plans x-critiques et des plans x-critiques "en apparence".
## gamm est le tableau des polynômes dont les racines sont des vraies singularités
## betaa est le tableau des paramétrisations rationnelles associées au tableau gamm
## lamda est le tableau des paramétrisations rationnelles associées au tableau gamm 
## AppSing est le tableau des polynômes dont les racines sont des singularités apparentes
## Appbeta est le tableau des paramétrisations rationnelles associées au tableau AppSing
## ATTENTION ses listes sont ordonnées
## F est la partie sans facteurs carrées du resultant
## SR est le tableau des coefficients sous-résultants
## L est la liste des polynomes sous-résultants donnée par la fonction StHa
## gam est la liste des polynomes gammas
## retourne une liste de deux listes : la première liste est la liste des valeurs 
## critiques et la deuxième est la liste des deuxième et troisième coordonnées correspondantes dans le 
## même ordre

CriticalFibers := proc(gamm,betaa,AppSing,Appbeta,L,gam,lamda,Applamda,F,SR,P,Q)
  local pc,i,j,co,n,k,m,PC,CO,Apc,Aco,APC,ACO,fibs,fib,h,f,Jac;

	
	n:=nops(L);
	k:=nops(gam);
	pc:=[];
	co:=[];
	fibs := [];
	Jac:=op(NullSpace(Jacobian([P,Q],[x,y,z])));
 for i from 1 by 1 to k do 
    for j from 1 by 1 to n  do
     if degree(gamm[i][j],x)>0 then
	PC[i][j] := [fsolve(gamm[i][j],x)];
	pc:=[op(pc),op(PC[i][j])];
	CO[i][j] :=[seq([eval(betaa[i][j],x=PC[i][j][m]),[eval(lamda[i][j],x=PC[i][j][m])]],m=1..nops(PC[i][j]))];
	co:=[op(co),op(CO[i][j])];
	f := deflate(F,gamm[i][j],betaa[i][j],i);
	fib:=[];
      for h from 1 to nops(PC[i][j]) do
	fib := [op(fib),LiftCritFib(CritFib(f,PC[i][j][h],CO[i][j][h][1],CO[i][j][h][2]),PC[i][j][h],SR)];
      od;
	fibs := [op(fibs),op(fib)];
     fi;
    od;
 od;

	Apc:=[];
	Aco:=[];

 for i from 1 by 1 to k do 
    for j from 2 by 1 to n  do
     if degree(AppSing[i][j],x)>0 then
	APC[i][j] := [fsolve(AppSing[i][j],x)];
	Apc:=[op(Apc),op(APC[i][j])];
	ACO[i][j] :=[seq(AppSingTri(APC[i][j][m],Jac,[eval(Appbeta[i][j],x=APC[i][j][m]),[fsolve(eval(Applamda[i][j],x=APC[i][j][m]),z)]]),m=1..nops(APC[i][j]))];
	Aco:=[op(Aco),op(ACO[i][j])];
	f := deflate(F,AppSing[i][j],Appbeta[i][j],i);
	fib:=[];
      for h from 1 to nops(APC[i][j]) do
	fib := [op(fib),LiftCritFib(CritFib(f,APC[i][j][h],ACO[i][j][h][1],ACO[i][j][h][2]),APC[i][j][h],SR)];
      od;
	fibs := [op(fibs),op(fib)];
     fi;
    od;
 od;

	pc:=[op(pc),op(Apc)];
 RETURN([pc,fibs]);
end:



## Calcul des deuxième et troisième coordonnées des points d'un plan régulier
## SR est le tableau des coefficients sous-résultants
## del est la liste des polynomes deltas
## v est la valeure régulière
ComputeRegFiber := proc(SR,del,v)
 local i,j, Y,Z,YY,ZZ,f;
	
	Y:=[];
	Z:=[];

 for i from 1 by 1 to nops(del) do
   if degree(del[i]) > 0 then
	
	YY:=[];
	YY:=[fsolve(subs(x=v,del[i]),y)];
	f:=-SR[i][i-1]/(i*SR[i][i]);
	f:=eval(f,x=v);
	ZZ:=[seq(eval(f,y=YY[j]),j=1..nops(YY))];
 	Y:=[op(Y),op(YY)];
	Z:=[op(Z),op(ZZ)];
   fi;
 od;
RETURN(MonTri([Y,Z]));
end:


## connection d'un plan critique à un plan  régulier
LRConnectCritToReg := proc(C,R)
  local Up, reg, Down, d, r, i, j, u,yUp,zUp,yDown,zDown,yreg,zreg,k;

	Up:=C[2];
	yUp:=Up[1];
	zUp:=Up[2];
	u:=nops(yUp);
	Down:=C[4];
	yDown:=Down[1];
	zDown:=Down[2];
	d := nops(yDown);
	reg := R[2];
	yreg:=reg[1];
	zreg:=reg[2];
 	r := nops(yreg);
	j := r - d - u;
	k:=nops(C[3][2]);


  if (k=1) then

  RETURN([seq([ [C[1],yUp[i],zUp[i]],[R[1],yreg[i],zreg[i]] ],i=1..u),seq([ [C[1],C[3][1],C[3][2][1]],[R[1],yreg[u+i],zreg[u+i]] ],i=1..j), seq([ [C[1],yDown[d-i+1],zDown[d-i+1]],[R[1],yreg[r-i+1],zreg[r-i+1]] ],i=1..d)]);
  else

 RETURN([seq([ [C[1],yUp[i],zUp[i]],[R[1],yreg[i],zreg[i]] ],i=1..u),seq([ [C[1],C[3][1],C[3][2][k-i+1]],[R[1],yreg[u+i],zreg[u+i]] ],i=1..k), seq([ [C[1],yDown[d-i+1],zDown[d-i+1]],[R[1],yreg[r-i+1],zreg[r-i+1]] ],i=1..d)])
fi;

end:






## connection d'un plan critique à un plan  régulier
RLConnectCritToReg := proc(C,R)
  local Up, reg, Down, d, r, i, j, u,yUp,zUp,yDown,zDown,yreg,zreg,k;

	Up:=C[2];
	yUp:=Up[1];
	zUp:=Up[2];
	u:=nops(yUp);
	Down:=C[4];
	yDown:=Down[1];
	zDown:=Down[2];
	d := nops(yDown);
	reg := R[2];
	yreg:=reg[1];
	zreg:=reg[2];
 	r := nops(yreg);
	j := r - d - u;
	k:=nops(C[3][2]);

  if (k=1) then

  RETURN([seq([ [C[1],yUp[i],zUp[i]],[R[1],yreg[i],zreg[i]] ],i=1..u),seq([ [C[1],C[3][1],C[3][2][1]],[R[1],yreg[u+i],zreg[u+i]] ],i=1..j), seq([ [C[1],yDown[d-i+1],zDown[d-i+1]],[R[1],yreg[r-i+1],zreg[r-i+1]] ],i=1..d)]);
  else

 RETURN([seq([ [C[1],yUp[i],zUp[i]],[R[1],yreg[i],zreg[i]] ],i=1..u),seq([ [C[1],C[3][1],C[3][2][i]],[R[1],yreg[u+i],zreg[u+i]] ],i=1..k), seq([ [C[1],yDown[d-i+1],zDown[d-i+1]],[R[1],yreg[r-i+1],zreg[r-i+1]] ],i=1..d)])
fi;

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:



GenericTest:=proc(V,W)
 local D,SR,del,F,L,cond,G,gam,sr,beta,IG,AppSing,Appbeta,Applamda,gamm,betaa,lamda,p,q,s,g,P,Q;

	P:=expand(V);
	Q:=expand(W);
	p:=degree(P,z);
    if (p > 0) then 
	s:=lcoeff(P,[z]);
	cond:=type(s,realcons);
    else
	cond:=true;
    fi;

    if (cond = true) then
	q:=degree(Q,z);
if (q > 0) then
	g:=lcoeff(Q,[z]);
	cond:=type(g,realcons);
    fi;
    fi;

    if (cond = false) then
	RETURN(cond);
    else
	D:=Delt(P,Q);
	SR:=eval(D[1]);
	del:=D[2];
	F:=D[3];
	L:=D[4];
	cond:=IsPseudoGeneric(SR,del);
    fi;
 
  if (cond = false) then
	RETURN(cond);
  else 
	G := Gam(F);
	gam := G[2];
	sr := G[1];
	cond:=IsGeneric(sr,gam);
   if (cond = false) then 
	RETURN(cond)
   else
	beta := ComputeBetaCrit(sr,gam);
	IG:=IsGGeneric(SR,gam,beta,L);
	cond:=IG[1];
    if (cond = false) then
	RETURN(cond);
    else
	AppSing:=IG[2];
	Appbeta:=IG[3];
	Applamda:=IG[4];
	gamm:=IG[5];
	betaa:=IG[6];
	lamda:=IG[7];
	RETURN(cond);
    fi;
   fi;
  fi;
end:
		




MaTopo := proc(P,Q)
  local IG, D, G,gam, sr, cf, beta, SR, del,F, L, cv, complexe, rv, rf,C,R,AppSing,Appbeta,Applamda,gamm,betaa,lamda;

  #calcul des fibre critiques
	
	D:=Delt(P,Q);
	SR:=eval(D[1]);
	del:=D[2];
	F:=D[3];
	L:=D[4];


	G := Gam(F);
	gam := G[2];
	sr := G[1];
	beta := ComputeBetaCrit(sr,gam);

	IG:=IsGGeneric(SR,gam,beta,L);
	AppSing:=IG[2];
	Appbeta:=IG[3];
	Applamda:=IG[4];
	gamm:=IG[5];
	betaa:=IG[6];
	lamda:=IG[7];

	cf:=CriticalFibers(gamm,betaa,AppSing,Appbeta,L,gam,lamda,Applamda,F,SR,P,Q);
	cf:=MonTri(cf);
	cv := cf[1];
	cf := cf[2];
 
 ## calcul des fibres régulières

	rv := RegValues(cv);
	rf := [seq(ComputeRegFiber(SR,del,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(LRConnectCritToReg(C[i],R[i+1]),i=1..nops(C)),seq(RLConnectCritToReg(C[nops(C)-i+1],R[nops(C)-i+1]), i=1..nops(C))];
  RETURN(complexe); 
end:

Draw := proc(P,Q)
  local comp,C,i;

	comp := MaTopo(P,Q);
	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: