#Function showpaths

#Just a function which plot the minimum spanning tree of the critical
# points of F.

showpaths:=proc(F,incx,incy)
local crit,tree,pol,x,y,RF;
    crit:=table();
    tree:=table();
    pol:=table();
    x,y:=incx,incy;
    RF:=sqrfree(resultant(F,diff(F,y),y))[2];
    `if`(nargs=3,pointstree(RF,x,crit,pol),pointstree(RF,x,crit,pol,args[4]));
    findtree(crit,tree);
#    intermediarypoints(crit,tree);
    plotgraph(crit,tree);
end;

#Function onlypatha

# Here we only compute the paths. This function is to test this part.

onlypaths:=proc(F,x,y)
local crit,tree,pol,t1,RF,sqrRF;
    t1:=time();
# we first define table used in our program
    crit:=table();
    pol:=table();
    tree:=table();
    RF:=evala(resultant(numer(F),diff(numer(F),y),y));
    sqrRF:=sqrfree(RF,x)[2];
    `if`(nargs=3,pointstree(sqrRF,x,crit,pol),pointstree(sqrRF,x,crit,pol,args[4]));
	findtree(crit,tree);
    localorder(crit);
    findpath(crit,tree);
    crit,tree,time()-t1;;
end;

# Function monodromy

# That's the general function.

# Input: - F the function which defines the Riemann surface.
#        - x the first unknown.
#        - y the second unknown.
# Output: - a the base point.
# 	      - ya the points of the Riemann surface over a.
#    	  - mon the monodromy group, which is presented as a list of couple b[i],permurtation group associated at the critical point.

monodromy:=proc(f,incx,incy)
local RF,F,x,y,t1,t0,i,numtot,con,monodro,alpha,k,T,a,basept,showtree,base,optio,opt,pol,isdebug,crit,d,poly,tree,sqrRF,z,K,M,SR,ext,j,prime,prov,S;
    t1:=time();
# we first define table used in our program
    crit:=table();
    pol:=table();
    tree:=table();
# now we look at options which are in the input.
    optio:=[args[4..nargs]]; opt:=map2(op,1,optio);
    if member('primenumber',opt,i) then
        prime:=op([i,2],optio)
    fi;
    if member('basepoint',opt,'i') then
        basept:=true;
        base:=op([i,2],optio);
    else
        basept:=false;
    fi;
    showtree:=member('displaygraph',opt);
#    if member('numstep',opt,'i') then
#        step:=true;
#        nstep:=op([i,2],optio);
#    else
#        step:=false;
#    fi;
    isdebug:=member('debugging',opt);
# we fix F with local unknowns.
    F:=eval(f,{incx=x,incy=y,I=RootOf(x^2+1)});
    d:=degree(F,y);
# we now begin computations
    t0:=time();
# We want only one primitive element for the coefficient field. Thus,
# we begin to create this primitive element if necessary.
    ext:=`algcurves/g_ext_r`(F);
    if nops(ext)>1 then
        K:=evala(Primifield(F));
        F:=eval(F,K[2]);
        M:=eval(op([1,1,1,1],K),_Z=z);
    elif nops(ext)=1 then
        M:=eval(op([1,1],ext),_Z=z);
        F:=evala(F);
    else
        M:=z;
    fi;
    RF:=evala(resultant(numer(F),diff(numer(F),y),y));
    sqrRF:=sqrfree(RF,x)[2];
    if not(assigned(prime)) then
        SR:=evala(RF/gcd(RF,diff(RF,x)));
        S:=evala(resultant(SR,diff(SR,x),x));
        prime:=LVGP(F,x,y,M,z,S);
    fi;
 userinfo(3,'monodromy',`Computing the points of the tree`);
    `if`(basept,pointstree(sqrRF,x,crit,pol,base),pointstree(sqrRF,x,crit,pol));
 userinfo(3,'monodromy',`running time`,time()-t0);
    t0:=time();
 userinfo(3,'monodromy',`Computing the tree`);
	findtree(crit,tree);
 userinfo(3,'monodromy',`running time`,time()-t0);
    t0:=time();
 userinfo(3,'monodromy',`Computing singular parts of Puiseux expansions`);
    poly:=pol['list'];
    for i in poly do
       prov:=mypuiseux(F,x,y,i[1],0,'valuation'=i[2],'primenumber'=prime,T,'roots'=pol['roots'][i[1]],'moreterms'
#,'irreducible'
);
        for j in prov do
            crit[pol[op([1,1,2,2],j)]['index']]['pui']:=j;
#print(j,"coucou",seq(op([1,2,2],j),k in j));
            crit[pol[op([1,1,2,2],j)]['index']]['ramif']:=[seq(degree(op([1,2],k),T),k in j)];
#print(pol[op([1,1,2,2],j)]['index']);
        od;
    od;
# return([crit,pol,tree,[alpha,T,x,y]]);
userinfo(3,'monodromy',`running time`,time()-t0);
    t0:=time();
 userinfo(3,'monodromy',`Computing intermediary points`);
    numtot:=intermediarypoints(crit,tree);
 userinfo(3,'monodromy',`running time`,time()-t0);
 userinfo(2,'monodromy',`Total number of intermediary points`,numtot);
	if showtree then plotgraph2(crit,tree); fi;
    t0:=time();
 userinfo(3,'monodromy',`Computing the precision needed for evaluations`);
    precision(F,x,y,d,tree,crit['nbcrit']);
 userinfo(3,'monodromy',`running time`,time()-t0);
    t0:=time();
 userinfo(3,'monodromy',`Computing truncation orders`);
    truncationorder(F,x,y,d,crit,tree);
 userinfo(3,'monodromy',`running time`,time()-t0);
    t0:=time();
 userinfo(3,'monodromy',`Computing Puiseux expansions`);
    setpuiseux(F,x,y,alpha,T,crit,pol,tree);
 userinfo(3,'monodromy',`running time`,time()-t0);
    t0:=time();
 userinfo(3,'monodromy',`Computing local orders`);
#return(crit);
    localorder(crit);
 userinfo(3,'monodromy',`running time`,time()-t0);
    t0:=time();
 userinfo(3,'monodromy',`Computing the set of paths`);
    findpath(crit,tree);
 userinfo(3,'monodromy',`running time`,time()-t0);
    t0:=time();
 #userinfo(3,'monodromy',`Computing the list of connection we have to do for each edge of the tree`);
 #   segconnect(crit,tree);
 #userinfo(3,'monodromy',`running time`,time()-t0);
 #   t0:=time();
 userinfo(3,'monodromy',`Computing the connection of each vertex`);
    con:=link(crit,tree,alpha,x,d);
 userinfo(3,'monodromy',`running time`,time()-t0);
    t0:=time();
 userinfo(3,'monodromy',`Computing the monodromy`);
    monodro:=mono(con,crit,tree,d,x,alpha);
 userinfo(3,'monodromy',`running time`,time()-t0);
 userinfo(1,'monodromy',`total time`,time()-t1);
    a:=crit[0]['val'];
#    res:=[a,eval(crit[0]['pui'],x=a),[seq([crit[i]['val'],convert(monodro[i],'disjcyc')],i=1..crit['nbcrit'])]]
    `if`(isdebug,[
[a,eval(crit[0]['pui'],x=a),[seq([crit[i]['val'],convert(monodro[i],'disjcyc')],i=1..crit['nbcrit'])]],
crit,pol,tree,con,monodro,[alpha,T,x,y]],[a,eval(crit[0]['pui'],x=a),[seq([crit[i]['val'],convert(monodro[i],'disjcyc')],i=1..crit['nbcrit'])]])
#    `if`(isdebug,[res,[crit,pol,tree]],res);
end:

# Function esqrt

# Input: - z a complex number
#        - e the exponent
#        - alpha the argument of the branch cut
# Output: - the e-th root of z according to the branch cut defined by alpha

esqrt:=proc(z,e,alpha)
local evaz;
    evaz:=evalf(z);
    if evaz::complex and alpha::realcons then
        if evaz=0 then return 0 fi;
        if argument(evaz)<evalf(alpha) then
            evalf(abs(evaz)^(1/e)*exp(I*(argument(evaz)+2*Pi)/e))
        else
            evalf(abs(evaz)^(1/e)*exp(I*(argument(evaz))/e))
        fi;
    else
        esqrt(args):='esqrt'(args);
    fi;
end:

# Function pointstree

# This function computes the set of the critical points and the point base, and then sorts the critical points.
# Input: - Delta the discriminant of F multiplyed by the leading coefficient of F.
#        - x the unknown of Delta.
#        - data the module
#   	 - c (optionnal) the base point for the monodromy.
# Output: all what we need is exported in the module data.


basepoint:=proc(b)
local i,num,real,dist;
    real:=Re(b[1]);
    num:=1;
    for i from 2 to nops(b) do
        if Re(b[i])<real then
            real:=Re(b[i]);
            num:=i;
        fi;
    od;
    dist:=`if`(num<>1,abs(b[1]-b[num]),abs(b[1]-b[2]));
    for i from 2 to nops(b) do
        dist:=`if`(i<>num and abs(b[i]-b[num])<dist,abs(b[i]-b[num]),dist);
    od;
    b[num]-dist;
end:

# basepoint crée un point de base "intelligent" par rapport à ma
# stratégie de nombre d'étapes.
# Problème: il ne serait pas mieux que le point de base soit réel ?

pointstree:=proc(RF,x,crit,pol)
local p,numpoints,prov,b,a,ord,i,j,k;
    p:=nops(RF);
    pol['list']:=RF;
# rootof is the set of critical points not indexed.
    numpoints:=NULL;
    for i in RF do
        prov:=[fsolve(i[1]=0,x,'complex')];
        pol['roots'][i[1]]:=prov;
        pol['dig'][i[1]]:=Digits;
        numpoints:=numpoints,prov;
    od;
    numpoints:=[numpoints];
    prov:=table();
    for i to p do
        for j in numpoints[i] do
            prov[j]:=RF[i,1];
        od;
    od;
    b:=[seq(op(numpoints[i]),i=1..p)];
    crit['nbcrit']:=nops(b);
# b is the set of critical points expressed in numeric.
    a:=`if`(nargs=5,evalf(args[5]),basepoint(b));
    crit[0]['val']:=a;
# a is the base point.
    ord:=proc(b1,b2) evalb((argument(b1-a)<argument(b2-a)) or ((argument(b1-a)=argument(b2-a)) and (abs(b1-a)<abs(b2-a)))); end;
    b:=sort(b,ord);
# we sort critical points.
    crit['list']:=b;
    for i to nops(b) do
        crit[i]['dig']:=Digits;
        crit[i]['val']:=b[i];
        crit[i]['pol']:=prov[b[i]];
        pol[b[i]]['index']:=i;
    od;
    for i to p do
        pol[RF[i,1]]['index']:=[seq(`if`(member(k,b,'j'),j,NULL),k in numpoints[i])];
    od;
# it links the factors of the discriminant and their numerical roots.
    i:=1;
    while Im(b[i])<Im(a) do i:=i+1 od;
    crit['limit']:=i;
end:


# Function findtree
# It computes the minimal spanning tree which links the critical points and the base point.

findtree:=proc(crit,tree)
    local pts,seen,notseen,numtree,value,arg,numseen,numnotseen,i,j,better,numbseen,numbnotseen;
    pts:=[seq(eval(crit[i]['val']),i=1..crit['nbcrit'])];
	seen:={crit[0]['val']};
	notseen:={op(pts)};
    numtree:=1:
	while notseen<>{} do
		value:=abs(seen[1]-notseen[1]):
		arg:=argument(seen[1]-notseen[1]):
		numseen:=1:
		numnotseen:=1:
		for i to nops(notseen) do
			for j to nops(seen) do
                better:=abs(seen[j]-notseen[i])<value
                or (abs(seen[j]-notseen[i])=value and  argument(seen[j]-notseen[i])<arg)
                or (abs(seen[j]-notseen[i])=value and  argument(seen[j]-notseen[i])=arg)
                    and (Re(notseen[i])<Re(notseen[numnotseen])
                         or (Re(notseen[i])=Re(notseen[numnotseen]) and Im(notseen[i])<Im(notseen[numnotseen])));
# we search a point nearest from the tree than notseen[numnotseen]
                if better then
                    value:=abs(seen[j]-notseen[i]);
					arg:=argument(seen[j]-notseen[i]);
					numseen:=j;
					numnotseen:=i;
				fi:
			od:
		od:
        numbseen:=`if`(member(seen[numseen],pts,'i'),i,0);
        member(notseen[numnotseen],pts,'numbnotseen');
        tree[numtree]['seg']:=[numbseen,numbnotseen];
        crit[numbnotseen]['found']:=numtree;
        tree[[numbseen,numbnotseen]]['index']:=numtree;
        tree[[numbnotseen,numbseen]]['index']:=numtree;
        if assigned(crit[numbseen]['close']) then
            crit[numbseen]['close']:=crit[numbseen]['close'],numbnotseen;
        else
            crit[numbseen]['close']:=numbnotseen;
        fi;
        crit[numbnotseen]['close']:=numbseen;
		numtree:=numtree+1;
		seen:= seen union {pts[numbnotseen]};
		notseen:=notseen minus {pts[numbnotseen]};
	od;
    for i from 0 to crit['nbcrit'] do
        crit[i]['close']:=[crit[i]['close']];
    od;
end;


# Function intermediarypoints

# This function computes all intermediary points needed to compute the
# monodromy (points where we have to compute puiseux expansions and
# points where we evaluate them).

intermediarypoints:=proc(crit,tree)
local numtot,step,i,segm,pt1,pt2,mi,delta,dm,p,j,prov,nextseg,nb;
    step:=proc() ceil(evalf(ln(dm/delta)/ln(2))) end;
# we look the segments of the tree one by one.
    nb:=crit['nbcrit'];
    numtot:=nb:
# numtot is the number of numeric puiseux expansions we will have to
# do: we first have all midpoints.
    crit[0]['nearcrit']:=crit[tree[1]['seg'][2]]['val'];
	for i to nb do
		segm:=tree[i]['seg'];
        pt1:=crit[segm[1]]['val'];
        pt2:=crit[segm[2]]['val'];
        mi:=1/2*(pt1+pt2);
		tree[i]['mid']:=mi;
# we first look the segment [pt1,m].
        delta:=abs(pt1-crit[segm[1]]['nearcrit']);
		dm:=abs(pt1-mi);
        p:=step();
        prov:=[seq(2^j*delta/dm*mi+(1-2^j*delta/dm)*crit[segm[1]]['val'],j=0..p-1)];
        tree[i]['m'][1]:=prov;
        tree[i]['ev'][1]:=`if`(p<=0,[(pt1+mi)/2],[(pt1+prov[1])/2,seq((prov[j]+prov[j+1])/2,j=1..p-1),(prov[p]+mi)/2]);
# the intermediary points in m1[i] are sorted: the first is the nearest intermediary point from pt1 and the last is the further one.
        numtot:=numtot+max(0,p);
# we now look the segment [m[i],pts[segm[2]]].
        nextseg:=tree[i+1]['seg'];
        crit[segm[2]]['nearcrit']:=`if`(i=nb or nextseg[1]<>segm[2] or abs(pt1-pt2)<abs(crit[nextseg[2]]['val']-pt2),crit[segm[1]]['val'],crit[nextseg[2]]['val']);
        delta:=abs(pt2-crit[segm[2]]['nearcrit']);
# considere le point de base comme un point critique. Faudrait faire
# en sorte que non... a voir...
		dm:=abs(pt2-mi);
        p:=step();
        prov:=[seq(2^j*delta/dm*mi+(1-2^j*delta/dm)*crit[segm[2]]['val'],j=0..p-1)];
        tree[i]['m'][2]:=prov;
        tree[i]['ev'][2]:=`if`(p<=0,[(pt2+mi)/2],[(pt2+prov[1])/2,seq((prov[j]+prov[j+1])/2,j=1..p-1),(prov[p]+mi)/2]);
# the intermediary points in m2[i] are sorted: the first is the nearest intermediary point from pt2 and the last is the further one.
        numtot:=numtot+max(0,p);
    od;
    numtot;
end:

# Function plot_graph2
# It creates the graph of the intermediary points

plotgraph2:=proc(crit,tree)
local n,graph,P,pts,P1,P2,i;
    n:=crit['nbcrit'];
    pts:=seq([crit[tree[i]['seg'][1]]['val'],crit[tree[i]['seg'][2]]['val']],i=1..n);
    graph:=[seq(map(x->[Re(x),Im(x)],pts[i]),i=1..n)];
	P:=plot(graph,'axes'='BOXED','title'="Expansion and evaluation points along the euclidean minimum spanning tree");
    pts:=[seq(op(tree[i]['ev'][1]),i=1..n),seq(op(tree[i]['ev'][2]),i=1..n)];
    graph:=[seq([Re(pts[i]),Im(pts[i])],i=1..nops(pts))];
    P1:=plot(graph,'axes'='BOXED','style'='point','color'='red');
    pts:=[seq(op(tree[i]['m'][1]),i=1..n),seq(op(tree[i]['m'][2]),i=1..n),seq(tree[i]['mid'],i=1..n),seq(crit[i]['val'],i=1..n)];
    graph:=[seq([Re(pts[i]),Im(pts[i])],i=1..nops(pts))];
    P2:=plot(graph,'axes'='BOXED','style'='point','color'='blue');
    print(plots['display'](P,P1,P2));
end:

plotgraph:=proc(crit,tree)
#It creates the graph of the tree of critical points.
local n,pts,graph,i;
    n:=crit['nbcrit'];
    pts:=seq([crit[tree[i]['seg'][1]]['val'],crit[tree[i]['seg'][2]]['val']],i=1..n);
    graph:=[seq(map(x->[Re(x),Im(x)],pts[i]),i=1..n)];
    plot(graph,'axes'='BOXED');
end:


# Function precision
# It computes the set of precision needed in order to have reliable connections.
prec:=proc(F,x,y,d,x0,tree)
local provDig,inter,yy,rho,l,h,k;
    provDig:=Digits;
    inter:=true;
    while inter do
        yy:=fsolve(subs(x=x0,F)=0,y,complex);
        rho:=[seq(evalf(d*abs(subs({x=x0,y=yy[l]},F)/(mul(yy[l]-yy[k],k=1..l-1)*mul(yy[l]-yy[k],k=l+1..d)))),l=1..d)];
# rho is the list of the Smith's radius. We now look if the circle have an intersection
        for l to d-1 do
            for h from l+1 to d do
                inter:=`if`(rho[l]+rho[h]<abs(yy[l]-yy[h]),false,true);
                if inter then break; fi;
            od;
            if inter then break; fi;
        od;
        Digits:=`if`(inter,Digits+5,provDig);
    od;
    tree[x0]['fiber']:=[yy];
    tree[x0]['dig']:=Digits;
    min(seq(seq(abs(yy[l]-yy[h])-rho[l]-rho[h],h=l+1..d),l=1..d-1));
end;

precision:=proc(F,x,y,d,tree,n)
local i,j,eva;
    for i to n do
        eva:=tree[i]['ev'][1];
        tree[i]['epsilon'][1]:=[seq(prec(F,x,y,d,eva[j],tree),j=1..nops(eva))];
        eva:=tree[i]['ev'][2];
        tree[i]['epsilon'][2]:=[seq(prec(F,x,y,d,eva[j],tree),j=1..nops(eva))];
    od;
end:

# Function truncationorder
# It computes the set of the truncation orders needed in order to have reliable connections.

radius:=proc(delta,dist) if delta=infinity then infinity else (0*dist+1*delta)/1 fi; end:

maxcircle:=proc(pol,center,x,d,r)
# it computes an upper bound for roots of F(x,y) in D(center,r)
local prov,borne4,borne8,i,l;
    prov:=seq(add(abs(coeff(subs(x=center+x,pol[i]),x,l))*r^l,l=0..degree(pol[i],x)),i=1..d);
    borne4:=max(seq(evalf((d*prov[i])^(1/i)),i=1..d));
    borne8:=max(seq(evalf(prov[i]/binomial(d,i))^(1/i)*(2^(1/d)-1)^(-1),i=1..d));
    min(borne4,borne8);
end:

#maxcircle:=proc(pol,c,x,d,r)
# nouvelle version avec series pour gerer le cas ou F n'est pas
# unitaire... faut voir si c'est sur.
#local prov,i,n,borne4,borne8,j;
#    prov:=[seq(0,i=1..d)];
#    for i to d do
#        n:=max(degree(numer(pol[i]),x),degree(denom(pol[i]),x))+1;
#        prov[i]:=[op(series(pol[i],x=c,n))];
#        prov[i]:=add(abs(prov[i,2*j-1])*r^prov[i,2*j],j=1..nops(prov[i])/2-1);
#    od;
#   borne4:=max(seq(evalf((d*prov[i])^(1/i)),i=1..d));
#    borne8:=max(seq(evalf(prov[i]/binomial(d,i))^(1/i)*(2^(1/d)-1)^(-1),i=1..d));
#    min(borne4,borne8);
#end:


truncation:=proc(eps,M,r,dist) if r=infinity then 1 else max(1,ceil(evalf((ln(eps*(r-dist)/(4*M*r)))/ln(dist/r)))-1) fi; end:
# Attention: j'ai viré le ceil

truncrit:=proc(eps,e,pt,delta,dist,d,x,poly)
# it computes truncation order above branch points.
local r,M;
    r:=radius(delta^(1/e),dist^(1/e));
    M:=maxcircle(subs(x=x^e+pt,poly),0,x,d,r);
    truncation(eps,M,r,dist^(1/e))/e;
#    truncation(eps,M,r,dist^(1/e));
end:

truncationorder:=proc(F,x,y,d,crit,tree)
local ad,poly,i,segm,p,z,j,dist,delta,r,M,k,mm,prov;
    ad:=coeff(F,y,d);
    poly:=[seq(coeff(F,y,d-i)/ad,i=1..d)];
# for each computation of the truncation order: let x0 be the point
# where we will compute the Puiseux expansion, x1 the nearest one
# where we will evaluate this Puiseux expansion and b the nearest
# critical point from x0. Then the following words in the algorithm
# are :
# -> dist=abs(x0-x1)
# -> delta=abs(x0-b)
    for k to 2 do
        for i to crit['nbcrit'] do
            segm:=tree[i]['seg'];
            mm:=tree[i]['m'][k];
            p:=nops(mm);
            z:=abs(crit[segm[k]]['val']-crit[segm[k]]['nearcrit']);
            for j to p do
                dist:=`if`(j=1,z/2,z*2^(j-3));
                delta:=z*2^(j-1);
                r:=radius(delta,dist);
                M:=maxcircle(poly,mm[j],x,d,r);
                tree[i]['n'][k][j]:=truncation(tree[i]['epsilon'][k][j],M,r,dist);
                dist:=`if`(j=p,abs(tree[i]['mid']-mm[p])/2,z*2^(j-2));
                r:=radius(delta,dist);
                M:=maxcircle(poly,mm[j],x,d,r);
                tree[i]['n'][k][j]:=max(truncation(tree[i]['epsilon'][k][j+1],M,r,dist),tree[i]['n'][k][j]);
            od;
# we look the midpoint m[i]
            delta:=abs(tree[i]['mid']-crit[segm[k]]['val']);
            dist:=`if`(p=0,delta/2,dist);
            r:=radius(delta,dist);
            M:=maxcircle(poly,tree[i]['mid'],x,d,r);
            prov:=truncation(tree[i]['epsilon'][k][p+1],M,r,dist);
            tree[i]['n']['m']:=`if`(k=1,prov,max(prov,tree[i]['n']['m']));
# we look the critical point pts[segm[k]]
            delta:=abs(crit[segm[k]]['val']-crit[segm[k]]['nearcrit']);
            dist:=`if`(p>=1,z/2,dist);
#print(crit[segm[k]]['ramif']);
            for j in {`if`(segm[k]=0,1,op(crit[segm[k]]['ramif']))} do
                prov:=truncrit(tree[i]['epsilon'][k][1],j,crit[segm[k]]['val'],delta,dist,d,x,poly);
                crit[segm[k]]['n']:=`if`(type(crit[segm[k]]['n'],realcons),max(crit[segm[k]]['n'],prov),prov);
            od;
        od;
    od;
end:


# Function setpuiseux
# It computes the set of the puiseux expansions needed to compute the monodromy.


# à retravailler en utilisant des développements
# symboliques/numériques pour les points ramifiés.

myeval:=proc(pui,T,x,alpha)
# it values a Puiseux expansion expressed with T in x, computing all
# the conjugates.
# input: a puiseux expansion in the form [x=a*T^e+b,y=g(T)]
# output: g(1/a*(x-b)^(1/e)) e times (one for each root of x^e-1).
local prov,e,oneroot,k;
    prov:=convert(op([2,2],pui),'horner',T);
    e:=degree(op([1,2],pui),T);
    oneroot:=evalf([seq(RootOf(_Z^e-1,index=k),k=1..e)]);
# je mets ca plutot qu'un fsolve pour etre sur de l'ordre.
    seq(eval(prov,T=oneroot[k]*esqrt((x-coeff(op([1,2],pui),T,0))/coeff(op([1,2],pui),T,e),e,alpha)),k=1..e);
end:


setpuiseux:=proc(F,x,y,alpha,T,crit,pol,tree)
local i,t0,puis,tab,j,l,k,p,poly,m1,m2,d,S,G,N,co,ind,pu;
    poly:=pol['list'];
    for i in poly do
        t0:=time():
        for j in pol[i[1]]['index'] do
            pu:=crit[j]['pui'];
            puis:=NULL;
            for k in pu do
                d:=degree(op([2,2],k),T);
                G:=evalf(expand(eval(F,{x=eval(op([1,2],k),T=x),y=eval(op([2,2],k),T=x)+x^d*y})));
                co:=[coeffs(G,{x,y},'ind')];
                ind:=[ind];
                G:=add(`if`(degree(ind[l],x)<k[3],0,co[l]*ind[l]/x^(k[3])),l=1..nops(ind));
                N:=crit[j]['n'];
                S:=regular(G,diff(G,y),x,y,0,degree(op([1,2],k),T)*N-d,true);
                S:=[op(1,k),y=op([2,2],k)+T^d*(eval(convert(S,'horner'),x=T))];
                puis:=puis,S;
            od;
            puis:=[puis];
            tab:=[seq(degree(op([1,2],p),T),p in puis)];
            crit[j]['perm']:=[seq(op([seq(add(tab[k],k=1..p-1)+l+1,l=1..tab[p]-1),add(tab[k],k=1..p-1)+1]),p=1..nops(tab))];
            crit[j]['pui']:=[seq(myeval(l,T,x,alpha), l in puis)];
        od;
   od;
 userinfo(4,'monodromy',`running time`,time()-t0);
# now we compute numerical puiseux expansion above all others points.
 userinfo(4,'monodromy',`computing puiseux expansions over`,crit[0]['val']);
 userinfo(4,'monodromy',`the developement is done with truncation`,crit[0]['n']);
    t0:=time();
    crit[0]['pui']:=allseries(F,x,y,crit[0]['val'],crit[0]['n']+1):
 userinfo(4,'monodromy',`running time`,time()-t0);
    for i to crit['nbcrit'] do
 userinfo(4,'monodromy',`computing expansions over`,tree[i]['mid']);
 userinfo(4,'monodromy',`the developement is done with truncation`,tree[i]['n']['m']);
        t0:=time();
        tree[i]['pui']['m']:=allseries(F,x,y,tree[i]['mid'],tree[i]['n']['m']+1):
 userinfo(4,'monodromy',`running time`,time()-t0);
        m1:=tree[i]['m'][1];
        for j to nops(m1) do
 userinfo(4,'monodromy',`computing expansions over`,m1[j]);
 userinfo(4,'monodromy',`the developement is done with truncation`,tree[i]['n'][1][j]);
            t0:=time();
            tree[i]['pui'][1][j]:=allseries(F,x,y,m1[j],tree[i]['n'][1][j]+1):
 userinfo(4,'monodromy',`running time`,time()-t0);
        od;
        m2:=tree[i]['m'][2];
        for j to nops(m2) do
 userinfo(4,'monodromy',`computing puiseux expansions over`,m2[j]);
 userinfo(4,'monodromy',`the developement is done with truncation`,tree[i]['n'][2][j]);
            t0:=time();
            tree[i]['pui'][2][j]:=allseries(F,x,y,m2[j],tree[i]['n'][2][j]+1):
 userinfo(4,'monodromy',`running time`,time()-t0);
        od;
    od;
 userinfo(4,'monodromy',`Done computing puiseux expansions`);
end:


allseries:=proc(F,x,y,pt,N)
local ptF,yy,i,dF;
    ptF:=numshift(F,x,pt);
    yy:=[fsolve(eval(ptF,x=0),y,'complex')];
    dF:=diff(ptF,y);
    [seq(eval(convert(regular(ptF,dF,x,y,i,N,true),'horner'),x=x-pt),i in yy)];
end;

###############################################
# NEW

localorder:=proc(crit)
local i,ord;
    ord:=proc(a,b) evalb(argument(crit[a]['val']-crit[i]['val'])<argument(crit[b]['val']-crit[i]['val'])) end;
    for i from 0 to crit['nbcrit'] do
        crit[i]['close']:=sort(crit[i]['close'],ord)
    od;
end;


tauplane:=proc(alpha,a,b)
# It answers true if alpha belongs to tau-plane defined by the edge
# Delta=[a,b] and false otherwise. The edge Delta belongs to a line
# which divides the plane in two parts. The tau-plane is defined as
# the half plane P such that alpha belongs to P if and only if
# arg(alpha-a)>=arg(b-a)
    if argument(b-a) <=0 then
        evalb(argument(alpha-a)>argument(b-a) and argument(alpha-a)<= evalf(Pi)+argument(b-a))
    else
        not(evalb(argument(alpha-a)<=argument(b-a) and argument(alpha-a)> evalf(-Pi)+argument(b-a)))
    fi;
end;

inter:=proc(tree,crit,l)
# It gives the edges of the tree which intersect the segment
# [a,alpha_l].
local L,i,a,b,f,t,T1,T2,A,B;
    f:=proc() (Im(crit[b]['val'])-Im(crit[a]['val']))*T1+(Re(crit[a]['val'])-Re(crit[b]['val']))*T2=Re(crit[a]['val'])*Im(crit[b]['val'])-Im(crit[a]['val'])*Re(crit[b]['val']) end;
    T1:=t*Re(crit[l]['val']-crit[0]['val'])+Re(crit[0]['val']);
    T2:=t*Im(crit[l]['val']-crit[0]['val'])+Im(crit[0]['val']);
    L:=NULL;
    A:=crit[0]['val'];
    B:=crit[l]['val'];
    for i to crit['nbcrit'] do
        a,b:=op(tree[i]['seg']);
        if not member(0,[a,b]) then
            if Re(B)>=Re(A) then
                if min(a,b) < l and l < max(a,b) and max(Re(crit[a]['val']),Re(crit[b]['val']))>max(Re(A)) then
                    L:=L,[i,fsolve(f(),t)];
                fi;
            elif Im(B)<Im(A) then
                if (min(a,b) < l and l < max(a,b) and Im(crit[max(a,b)]['val']) < Im(A)) or (min(a,b) > l and max(Re(crit[a]['val']),Re(crit[b]['val']))<Re(A) and Im(crit[max(a,b)]['val']) > Im(A) and Im(A) > Im(crit[min(a,b)]['val'])) then
                    L:=L,[i,fsolve(f(),t)];
                fi;
            else
                if (min(a,b) < l and l < max(a,b) and Im(crit[min(a,b)]['val']) > Im(A)) or (max(a,b) < l and max(Re(crit[a]['val']),Re(crit[b]['val']))<Re(A) and Im(crit[max(a,b)]['val']) > Im(A) and Im(A) > Im(crit[min(a,b)]['val'])) then
                    L:=L,[i,fsolve(f(),t)];
                fi;
            fi;
        fi;
    od;
# L contient les intersections avec la demi droite [a,alpha_l)
    L:=sort([L,['pt',1]],proc(i,j) evalb(i[2]<j[2]) end);
    i:=1;
    while i<=nops(L) and L[i,2]<=1 do i:=i+1 od;
    crit[l]['inter']:=L[1..i-1,1];
    crit[l]['moreinter']:=L[i..-1,1];
end;

connect:=proc(crit,tree,i,j,l)
# It gives the list of edges of the tree which connect the i-th edge
# to the j-th edge
# i may be 0 for the first step (it means that we start from the base
# point ; thus, we do not know which edge of the tree we have to take
# without considering the edge j.
local searched,path,k,path2;
    path:=NULL;
    searched:=`if`(j='pt',l,tree[max(i,j)]['seg'][1]);
    while searched <> 0 and `if`(i=0,true,not(member(searched,tree[`if`(j='pt',i,min(i,j))]['seg']))) do
        k:=crit[searched]['found'];
        path:=k,path;
        searched:=tree[k]['seg'][1];
    od;
    if searched=0 and i<>0 then
        searched:=tree[`if`(j='pt',i,min(i,j))]['seg'][1];
        path2:=NULL;
        while searched <> 0 do
            k:=crit[searched]['found'];
            path2:=path2,k;
            searched:=tree[k]['seg'][1];
        od;
        path:=path2,path;
    fi;
    [path];
end;


sens:=proc(tau,L)
    local eps,i;
    eps:=-1;
    for i in tau do
        if member(i,L) then
            eps:=-eps;
        fi;
    od;
    eps;
end;


followtree:=proc(crit,tree,i1,i2,h,l,eps)
local j,L,Lg,path,b,k,w,sigma,n,i,ww,a1,a2,a3;
    if h=0 then
        j:=0;
        L:=crit[0]['close'];
        if eps=1 then
            Lg:=seq(`if`(i>l,i,NULL),i in L);
            k:=`if`(Lg=NULL,min(op(L)),min(Lg));
        else
            Lg:=seq(`if`(i<l,i,NULL),i in L);
            k:=`if`(Lg=NULL,max(op(L)),max(Lg));
        fi;
        path:=[j,k];
    else
        if eps=-1 then
            if tauplane(crit[tree[i1]['seg'][1]]['val'],crit[0]['val'],crit[l]['val']) then
                j,k:=op(tree[i1]['seg']);
            else
                k,j:=op(tree[i1]['seg']);
            fi;
        else
            if tauplane(crit[tree[i1]['seg'][1]]['val'],crit[0]['val'],crit[l]['val']) then
                k,j:=op(tree[i1]['seg']);
            else
                j,k:=op(tree[i1]['seg']);
            fi;
        fi;
        path:=NULL;
    fi;
#    path:=i1;
    b:=true;
    while b do
        sigma:=crit[k]['close'];
        member(j,sigma,'n');
        if eps=1 then
            w:=`if`(n=nops(sigma),sigma[1],sigma[n+1])
        else
            w:=`if`(n=1,sigma[-1],sigma[n-1])
        fi;
        path:=path,[j,w,eps];
        if assigned(crit[k]['con']) then
            crit[k]['con']:=crit[k]['con'] union {[j,w,eps]};
        else
            crit[k]['con']:={[j,w,eps]};
        fi;
        if (i2='pt' and w=l) then
            sigma:=crit[w]['close'];
            member(k,sigma,'n');
            if eps=1 then
                ww:=`if`(n=nops(sigma),sigma[1],sigma[n+1])
            else
                ww:=`if`(n=1,sigma[-1],sigma[n-1])
            fi;
            a1:=argument(crit[k]['val']-crit[l]['val']);
            a2:=argument(crit[ww]['val']-crit[l]['val']);
            a3:=argument(crit[0]['val']-crit[l]['val']);
            if ww=k or (eps = -1 and `if`(a1<a3,a1<a2 and a2<a3,a2>a3 or a2<a1)) or (eps=1 and `if`(a1>a3,a1>a2 and a2>a3,a2<a1 or a2>a3)) then
                b:=false;
                path:=path,[k,w]
            else
                path:=path,[k,w];
                j,k:=k,w;
            fi;
        elif (tree[[k,w]]['index']=i2 and ((eps=-1 and tauplane(crit[w]['val'],crit[0]['val'],crit[l]['val'])) or (eps=1 and tauplane(crit[k]['val'],crit[0]['val'],crit[l]['val'])))) then
            b:=false;
        else
            path:=path,[k,w];
            j,k:=k,w;
        fi;
    od;
    path;
end;


findpath:=proc(crit,tree)
local n,l,pathi,L,h,lastpt,pathh,eps,tau;
    n:=crit['nbcrit'];
    for l to n do
        if member(l,crit[0]['close']) then
            crit[l]['path']:=[tree[crit[l]['found']]['seg']];
        else
#        path:=NULL;
            inter(tree,crit,l);
            L:=crit[l]['inter'];
            tau:=connect(crit,tree,0,L[1],l);
            eps:=sens(tau,[op(L[2..-1]),op(crit[l]['moreinter'])]);
            lastpt:=`if`(member(l,tree[tau[-1]]['seg']),tree[tau[-1]]['seg'][1],tree[tau[-1]]['seg'][2]);
            if tauplane(crit[lastpt]['val'],crit[0]['val'],crit[l]['val']) then
                eps:=-eps;
            fi;
            if member(l,[seq(op(tree[i]['seg']),i in tau[1..-2])]) then eps:=-eps fi;
            pathi:=followtree(crit,tree,0,L[1],0,l,eps);
            for h to nops(L)-1 do
                tau:=connect(crit,tree,L[h],L[h+1],l);
                eps:=sens(tau,[op(L[h+2..-1]),op(crit[l]['moreinter'])]);
                lastpt:=`if`(tau=[],op({op(tree[L[h]]['seg'])} intersect {op(tree[L[h+1]]['seg'])}),`if`(L[h+1]='pt',tree[tau[-1]]['seg'][1],`if`(L[h]<L[h+1],tree[tau[-1]]['seg'][2],tree[tau[1]]['seg'][1])));
                if tauplane(crit[lastpt]['val'],crit[0]['val'],crit[l]['val']) then
                    eps:=-eps;
                fi;
                if tau <> [] and member(l,[seq(op(tree[i]['seg']),i in tau[1..-2])]) then eps:=-eps fi;
                pathh:=followtree(crit,tree,L[h],L[h+1],h,l,eps);
                if [pathi][-1][2]=[pathh][1][1] then
                    pathi:=pathi,pathh
                else
                    pathi:=pathi,[[pathh][1][1],[pathi][-1][2]],pathh
                fi;
            od;
            crit[l]['path']:=[pathi];
        fi;
    od;
end;



# Function link.

# It computes the links between the solutions of F(x,y) over b[tree[i,1]] and those over b[tree[i,2]] for each i.


numbranch:=proc(valj,eps,valy2,d,z)
# it computes the y-value of valy2 corresponding to the y-value valj
# with precision eps.
local res,k;
    for k to d do if abs(valj-valy2[k])<eps/2 then res:=k; break; else res:=z; fi: od;
    res
end:

lastpt:=proc(tree,i,k,d,x,z)
local mk,link,valy1,valy2,evm,prov,j,l;
    mk:=tree[i]['m'][k];
    if nops(mk)=0 then
        link:=table([seq(l,l=1..d)]);
    else
        evm:=tree[i]['ev'][k];
        valy1:=eval(tree[i]['pui']['m'],x=evm[nops(evm)]);
        valy2:=eval(tree[i]['pui'][k][nops(mk)],x=evm[nops(evm)]);
        prov:=table([seq(numbranch(valy1[l],tree[i]['epsilon'][k][nops(evm)],valy2,d,z),l=1..d)]);
        for j from nops(mk)-1 by -1 to 1 do
            valy1:=eval(tree[i]['pui'][k][j+1],x=evm[j+1]);
            valy2:=eval(tree[i]['pui'][k][j],x=evm[j+1]);
            prov:=table([seq(numbranch(valy1[prov[l]],tree[i]['epsilon'][k][j+1],valy2,d,z),l=1..d)]);
        od;
        link:=prov;
    fi;
    link;
end:


link:=proc(crit,tree,alpha,x,d)
local Dig,i,segm,link1,link2,valy1,valy2,j,prov,numfrom,numto,ptfrom,ptto,theta1,theta2,sens,theta,l,prov1,prov2,link,z,pt,cote1,cote2,eva1,eva2,provback1,provback2,puis1,puis2;
    Dig:=Digits;
	for i to crit['nbcrit'] do
        segm:=tree[i]['seg']:
# we first link the midpoint respectively to the nearest intermediary
# point point from b[segm[1]] and b[segm[2]].
        link1[i]:=lastpt(tree,i,1,d,x,z);
        link2[i]:=lastpt(tree,i,2,d,x,z);
# we now link the two fibers close to b[segm[1]] and b[segm[2]] of the
# i-th edge of the tree.
        eva1:=tree[i]['ev'][1][1];
#        tree[eva1]['fiber']:=fsolve(eval(F,x=eva1),y,'complex');
        eva2:=tree[i]['ev'][2][1];
#        tree[eva2]['fiber']:=fsolve(eval(F,x=eva2),y,'complex');
        puis1:=`if`(nops(tree[i]['m'][1])=0,tree[i]['pui']['m'],tree[i]['pui'][1][1]);
        puis2:=`if`(nops(tree[i]['m'][2])=0,tree[i]['pui']['m'],tree[i]['pui'][2][1]);
        valy1:=eval(puis1,x=eva1);
        prov1:=table([seq(numbranch(valy1[link1[i][l]],tree[i]['epsilon'][1][1],tree[eva1]['fiber'],d,z),l=1..d)]);
        if member(z,[seq(prov1[l],l=1..d)]) then error "increase Digits" fi;
        provback1:=table([seq(op([1,2,l,2],prov1)=op([1,2,l,1],prov1),l=1..d)]);
        valy2:=eval(puis2,x=eva2);
        prov2:=table([seq(numbranch(valy2[link2[i][l]],tree[i]['epsilon'][2][1],tree[eva2]['fiber'],d,z),l=1..d)]);
        if member(z,[seq(prov2[l],l=1..d)]) then error "increase Digits" fi;
        provback2:=table([seq(op([1,2,l,2],prov2)=op([1,2,l,1],prov2),l=1..d)]);
        link[[segm[1],segm[2]]]:=table([seq(prov2[provback1[l]],l=1..d)]);
        link[[segm[2],segm[1]]]:=table([seq(prov1[provback2[l]],l=1..d)]);
    od;
# we have considered all connections along edges of the tree. We now
# look loops around critical points.
    for i from 0 to crit['nbcrit'] do
        pt:=crit[i]['val'];
#print(crit[i]['con']);
        for j in `if`(assigned(crit[i]['con']),crit[i]['con'],[]) do
            ptfrom:=j[1];
# the point from which we come
            ptto:=j[2];
# the point where we're going
            numfrom:=tree[[ptfrom,i]]['index'];
# the tree's number from which we come
            numto:=tree[[ptto,i]]['index'];
# the tree's number where we're going
            theta1:=argument(crit[ptfrom]['val']-pt);
            theta2:=argument(crit[ptto]['val']-pt);
            sens:=j[3];
            theta:=(theta1+theta2)/2;
            if (theta1<theta2 and sens=1) or (theta1>theta2 and sens=-1) then
                if theta>0 then
                    theta:=theta-Pi;
                else
                    theta:=theta+Pi
                fi;
            fi;
# theta is the argument of the branch cut we want to use.
# in order to get this branch cut, we have to look at the coefficient a
# in x in the function esqrt and to evaluate this function in
# theta+argument(a).

# Correction : with the numerical newton puiseux algorithm, a is equal
# to 1. Thus, we do not need to consider this here.
            if tree[numfrom]['seg'][1]=i then

#                nfrom:=tree[numfrom]['seg'][2];
                cote1:=1;
            else
#                nfrom:=tree[numfrom]['seg'][1];
                cote1:=2;
            fi;
            if tree[numto]['seg'][1]=i then
#                nto:=tree[numto]['seg'][2];
                cote2:=1;
            else
#                nto:=tree[numto]['seg'][1];
                cote2:=2;
            fi;
#            Digits:=crit[i]['dig'];
            eva1:=tree[numfrom]['ev'][cote1][1];
            valy1:=eval(crit[i]['pui'],{x=eva1,alpha=theta});
            eva2:=tree[numto]['ev'][cote2][1];
            valy2:=eval(crit[i]['pui'],{x=eva2,alpha=theta});
            prov1:=table([seq(numbranch(tree[eva1]['fiber'][l],tree[numfrom]['epsilon'][cote1][1],valy1,d,z),l=1..d)]);
            if member(z,[seq(prov1[l],l=1..d)]) then error "increase Digits" fi;
            link[j]:=table([seq(numbranch(valy2[prov1[l]],tree[numto]['epsilon'][cote2][1],tree[eva2]['fiber'],d,z),l=1..d)]);
            if member(z,[seq(link[j][l],l=1..d)]) then error "increase Digits" fi;
            Digits:=Dig;
        od;
    od;
    link;
end:



# Function mono
# It is the final function. It connects the segments of the path one by one

mono:=proc(link,crit,tree,d,x,alpha)
local i,pathi,n,eva,valy1,valy2,prov,z,provback,monodro,j,con,conback,cote,loop,l;
    for i to crit['nbcrit'] do
        pathi:=crit[i]['path'];
        n:=tree[pathi[1]]['index'];
        eva:=tree[n]['ev'][1][1];
        valy1:=eval(crit[0]['pui'],x=eva);
        valy2:=tree[eva]['fiber'];
        prov:=table([seq(numbranch(valy2[l],tree[n]['epsilon'][1][1],valy1,d,z),l=1..d)]);
        if member(z,[seq(prov[l],l=1..d)]) then error "increase Digits" fi;
# we first connect the base point to the first fiber.
        for j in pathi do
            prov:=table([seq(link[j][prov[l]],l=1..d)]);
        od;
# prov contains connections for the path from the base point to the
# last evaluation point close to alpha_i.
        provback:=table([seq(op([1,2,l,2],prov)=op([1,2,l,1],prov),l=1..d)]);
        n:=tree[pathi[-1]]['index'];
        cote:=`if`(tree[n]['seg'][1]=i,1,2);
        eva:=tree[n]['ev'][cote][1];
        valy1:=eval(crit[i]['pui'],{x=eva,alpha=-Pi});
        valy2:=tree[eva]['fiber'];
        con:=table([seq(numbranch(valy2[l],tree[n]['epsilon'][cote][1],valy1,d,z),l=1..d)]);
        if member(z,[seq(con[l],l=1..d)]) then error "increase Digits" fi;
        conback:=table([seq(op([1,2,l,2],con)=op([1,2,l,1],con),l=1..d)]);
        loop:=[seq(conback[crit[i]['perm'][con[l]]],l=1..d)];
        monodro[i]:=[seq(provback[loop[prov[l]]],l=1..d)];
    od;
    monodro;
end: