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: