findpath:=proc(x0,t)
# it defines one of the path we want to follow
    [t*x0,x0*exp(I*Pi*t/2),(1-t)*x0*exp(I*Pi/2)]
end;

fiber:=proc(F,x,y,x0)
# it just computes the fiber of the plane algebraic curve defined by F
# at the point x0
    fsolve(eval(F,x=x0),y,'complex')
end;

compose:=proc(a,b)
# it compute b o a, where a and b are two compatible table
local c,i;
    c:=op(op(a));
# this variable seems to be needed to have the correct result
    table([seq(op(1,i)=b[op(2,i)],i=c)]);
end;

# initial_step=1/8

allturns:=proc(F,x,y)
# the main function: it computes the permutations generated by
# analytic continuation of 3 turns in the complex x turn.
local A,B,C,D,E,dF,dF2,t,step,i,T,res,prov;
    T:=time();
 userinfo(2,'continuation',`Computing derivatives`);
    step:=1/8;
    A:=diff(F,x);
    B:=diff(F,y);
#    C:=diff(A,x);
#    D:=diff(A,y);
#    E:=diff(B,y);
#    dF:=evalf(collect(normal(-A/B),y));
    dF:=evalf(-A/B);
#    dF2:=evalf(collect(normal((-C+2*dF*D+dF^2*E)/B),y));
#    dF2:=(-C+2*dF*D+dF^2*E)/B;
    dF2:=0;
 userinfo(2,'continuation',`running time`,time()-T);
    res,prov:=oneturn(F,dF,dF2,x,y,t,step,2*exp(I*Pi/6));
    prov:=oneturn(F,dF,dF2,x,y,t,step,2*exp(2*I*Pi/3),prov);
    res:=res,prov[1],oneturn(F,dF,dF2,x,y,t,step,2*exp(7*I*Pi/6),prov[2])[1];
 userinfo(1,'continuation',`Total running time`,time()-T);
    [res]
end;

allturns2:=proc(F,x,y)
# the main function: it computes the permutations generated by
# analytic continuation of 3 turns in the complex x turn.
local A,B,C,D,E,dF,dF2,t,step,i,T,res,prov;
    T:=time();
 userinfo(2,'continuation',`Computing derivatives`);
    step:=1/8;
    A:=diff(F,x);
    B:=diff(F,y);
    C:=diff(A,x);
    D:=diff(A,y);
    E:=diff(B,y);
#    dF:=evalf(collect(normal(-A/B),y));
    dF:=evalf(-A/B);
    dF2:=evalf(-(C+2*dF*D+dF^2*E)/B);
 userinfo(2,'continuation',`running time`,time()-T);
    res,prov:=oneturn(F,dF,dF2,x,y,t,step,2*exp(I*Pi/6));
    prov:=oneturn(F,dF,dF2,x,y,t,step,2*exp(2*I*Pi/3),prov);
    res:=res,prov[1],oneturn(F,dF,dF2,x,y,t,step,2*exp(7*I*Pi/6),prov[2])[1];
 userinfo(1,'continuation',`Total running time`,time()-T);
    res
end;

inverse:=proc(t)
# it computes the inverse of a table
    local l;
    table([seq(op([1,2,l,2],t)=op([1,2,l,1],t),l=1..nops(op(op(t))))]);
end;

oneturn:=proc(F,dF,dF2,x,y,t,step,x0)
# it computes the analytic continuation along the path beginning at x0
local paths,perm,i,T,prov;
 userinfo(2,'continuation',`Beginning one turn`);
    T:=time();
    paths:=findpath(x0,t);
    perm:=`if`(nargs=9,inverse(args[9]),continuation(args[1..5],paths[1],args[6..7]));
 if nargs=9 then userinfo(2,'continuation',`Analytic continuation for the path`,paths[1],`already computed`) fi;
    for i in paths[2..3] do
        prov:=continuation(args[1..5],i,args[6..7]);
        perm:=compose(perm,prov)
    od;
 userinfo(2,'continuation',`Running time for this turn`,time()-T);
    perm,prov
end;

continuation:=proc(F,dF,dF2,x,y,path,t,step)
# This process the analytic continuation along a given path
local t0,perm,T,T1;
 userinfo(2,'continuation',`Beginning annalytic continuation for the path`,path);
    T:=time();
    t0:=0;
 userinfo(4,'continuation',`Analytic continuation between`,evalf(eval(path,t=t0)),`and`,evalf(eval(path,t=t0+step)));
    T1:=time();
    perm:=link(args[1..8],t0);
 userinfo(4,'continuation',`Running time`,time()-T1);
    while t0<1 do
        t0:=t0+step;
 userinfo(4,'continuation',`Analytic continuation between`,evalf(eval(path,t=t0)),`and`,evalf(eval(path,t=t0+step)));
        T1:=time();
        perm:=compose(perm,link(args[1..8],t0));
 userinfo(4,'continuation',`Running time`,time()-T1);
    od;
 userinfo(2,'continuation',`Running time for this path`,time()-T);
    perm
end;

link:=proc(F,dF,dF2,x,y,path,t,step,t0)
# It computes one step of the analytic continuation
local pt,yy,newpt,newfiber,newyy,app,appyy,outy,close,i,m,j,mm,d,root,perm,a,b;
#    if step<1/1024 then error `the stepsize becomes too small` fi;
#    perm:=table();
    pt:=evalf(eval(path,t=t0));
    yy:=`if`(nargs=11 and args[10]='first',args[11],[fiber(F,x,y,pt)]);
    newpt:=evalf(eval(path,t=min(t0+step,1)));
    newfiber:=`if`(nargs=11 and args[10]='last',args[11],{fiber(F,x,y,newpt)});
    newyy:=newfiber;
    app:=(newpt-pt)^2*eval(dF2,x=pt)/2+(newpt-pt)*eval(dF,x=pt)+y;
    appyy:=[seq(eval(app,y=i),i=yy)];
    outy:=NULL;
    close:=false;
    for i in appyy do
        m:=infinity;
        for j in newyy do
            mm:=abs(i-j);
            if mm<m then
                d,m,root:=m,mm,j
            elif mm<d then
                d:=mm
            fi;
        od;
        outy:=outy,root;
        newyy:=newyy minus {root};
        if m>d/10 then close:=true fi;
    od;
    outy:=[outy];
    if close then
 userinfo(3,'continuation',`decreasing the step to`,step/2);
        a:=procname(args[1..7],step/2,t0,'first',yy);
        b:=procname(args[1..7],step/2,t0+step/2,'last',newfiber);
        perm:=compose(a,b)
    else
        perm:=table([seq(yy[i]=outy[i],i=1..nops(yy))]);
    fi;
    perm
end;


###################################################################

plotfibers:=proc(F,x,y)
# It will plot the successives fibers along one segment for the two
# methods using respectively first and second derivatives
local T,A,B,C,D,E,dF,dF2,tt,t,step,path,points,P1,P2,P3;
    T:=time();
    t:=`if`(nargs>=5,args[5],tt);
    step:=1/2;
    A:=diff(F,x);
    B:=diff(F,y);
#    C:=diff(A,x);
#    D:=diff(A,y);
#    E:=diff(B,y);
    dF:=evalf(-A/B);
 userinfo(9,'continuation',`first derivative`,dF);
#    dF2:=evalf(-(C+2*dF*D+dF^2*E)/B);
    dF2:=0;
 userinfo(9,'continuation',`second derivative`,dF2);
    path:=`if`(nargs>=4,args[4],2*t);
    points:=plotcontinuation(F,dF,dF2,x,y,path,t,step,`if` (nargs>=6,args[6],Re));
    P1:=seq(plot(i,color=red),i=points[1]);
# first derivative
#    P2:=seq(plot(i,color=blue),i=points[2]);
# second derivative
    P3:=plot(points[3],color=green,style=point);
# fiber
 userinfo(1,'continuation',`Total running time`,time()-T);
    plots[display](P1,P3);
end;

plotcontinuation:=proc(F,dF,dF2,x,y,path,t,step,func)
# This process the analytic continuation along a given path in two
# ways (one and two derivatives), and store all successive fibers
local t0,points,T,T1,y0,prov;
 userinfo(2,'continuation',`Beginning annalytic continuation for the path`,path);
    T:=time();
    t0:=0;
 userinfo(4,'continuation',`Analytic continuation between`,evalf(eval(path,t=t0)),`and`,evalf(eval(path,t=t0+step)));
    T1:=time();
    y0:=[fiber(F,x,y,eval(path,t=0))];
    points:=[[seq([[func(eval(path,t=0)),func(i)]],i in y0)]$2,{}];
    prov:=plotlink(args[1..8],t0,'all',func,y0);
    points:=[join(points[1],prov[1]),join(points[2],prov[2]),points[3] union prov[3],prov[4]];
 userinfo(4,'continuation',`Running time`,time()-T1);
    while t0<1 do
#print(t0);
        t0:=t0+step;
 userinfo(4,'continuation',`Analytic continuation between`,evalf(eval(path,t=t0)),`and`,evalf(eval(path,t=t0+step)));
        T1:=time();
        prov:=plotlink(args[1..8],t0,'all',func,points[4]);
        points:=[join(points[1],prov[1]),join(points[2],prov[2]),points[3] union prov[3],prov[4]];
 userinfo(4,'continuation',`Running time`,time()-T1);
    od;
 userinfo(2,'continuation',`Running time for this path`,time()-T);
    points
end;

plotlink:=proc(F,dF,dF2,x,y,path,t,step,t0,state,func,fib1)
# It computes one step of the analytic continuation in two
# ways (one and two derivatives), and store all successive fibers
#
# state is 'all',1 or 2, depending of if we are making analytic
# continuation for the two strategies, or only one of them
#
# func is what we use to represent y (like Re, Im or something else)
local pt,yy,newpt,newfiber,newyy,app1,app2,appyy1,appyy2,close1,close2,i,m,j,mm,d,root,points,a,b,points1,points2,points3,outy1,outy2,outfiber1,outfiber2,outfiber;
    pt:=evalf(eval(path,t=t0));
    yy:=fib1;
    newpt:=evalf(eval(path,t=min(t0+step,1)));
    newfiber:=`if`(nargs=13,args[13],{fiber(F,x,y,newpt)});
    newyy:=newfiber;
    close1,close2:=false$2;
    if state in {'all',1} then
        app1:=(newpt-pt)*eval(dF,x=pt)+y;
        appyy1:=[seq(eval(app1,y=i),i=yy)];
        outy1:=NULL;
        outfiber1:=NULL;
        for i in appyy1 do
            m:=infinity;
            for j in newyy do
                mm:=abs(i-j);
                if mm<m then
                    d,m,root:=m,mm,j
                elif mm<d then
                    d:=mm
                fi;
            od;
            outy1:=outy1,i;
            outfiber1:=outfiber1,root;
            newyy:=newyy minus {root};
            if m>d/10 then close1:=true; break; fi;
        od;
        outy1:=[outy1];
        outfiber1:=[outfiber1];
    fi;
    if state in {2,'all'} then
        app2:=(newpt-pt)^2*eval(dF2,x=pt)/2+(newpt-pt)*eval(dF,x=pt)+y;
        appyy2:=[seq(eval(app2,y=i),i=yy)];
        outy2:=NULL;
        outfiber2:=NULL;
        newyy:=newfiber;
        for i in appyy2 do
            m:=infinity;
            for j in newyy do
                mm:=abs(i-j);
                if mm<m then
                    d,m,root:=m,mm,j
                elif mm<d then
                    d:=mm
                fi;
            od;
            outy2:=outy2,i;
            outfiber2:=outfiber2,root;
            newyy:=newyy minus {root};
            if m>d/10 then close2:=true; break; fi;
        od;
        outy2:=[outy2];
        outfiber2:=[outfiber2];
    fi;
    if close1 or close2 then
        if close1 and close2 then
 userinfo(3,'continuation',`decreasing the step to`,step/2,`for the two ways`);
            a:=procname(args[1..7],step/2,t0,'all',func,yy);
            b:=procname(args[1..7],step/2,t0+step/2,'all',func,a[4],newfiber);
            points1:=join(a[1],b[1]);
            points2:=join(a[2],b[2]);
        elif close1 then
 userinfo(3,'continuation',`decreasing the step to`,step/2,`for the first derivative`);
            a:=procname(args[1..7],step/2,t0,1,func,yy);
            b:=procname(args[1..7],step/2,t0+step/2,1,func,a[4],newfiber);
            points1:=join(a[1],b[1]);
            points2:=[seq([[func(newpt),func(i)]],i in outy2)];
        else
 userinfo(3,'continuation',`decreasing the step to`,step/2,`for the second derivative`);
            a:=procname(args[1..7],step/2,t0,2,func,yy);
            b:=procname(args[1..7],step/2,t0+step/2,2,func,a[4],newfiber);
            points1:=[seq([[func(newpt),func(i)]],i in outy1)];
            points2:=join(a[2],b[2]);
        fi;
        points3:=a[3] union b[3];
        outfiber:=b[4];
   else
        if state='all' and outfiber1<>outfiber2 then error "bad connection" fi;
        points3:={seq([func(pt),func(i)],i in yy)};
        if state in {1,'all'} then
            points1:=[seq([[func(newpt),func(i)]],i in outy1)];
            outfiber:=outfiber1;
        fi;
        if state in {2,'all'} then
            points2:=[seq([[func(newpt),func(i)]],i in outy2)];
            outfiber:=outfiber2;
        fi;
    fi;
    [points1,points2,points3,outfiber];
end;

join:=proc(a,b)
# it is used to join the lists of the analytic process
local i;
    [seq([op(a[i]),op(b[i])],i=1..nops(a))]
end;


###################################################

ordre:=proc(perms)
# This labels the fiber and gives the permutations in perms with their labelling.
local l,prov,n,lab,lab2,i,j;
    l:=nops(perms);
    prov:=op(op(perms[1]));
    n:=nops(prov);
    lab:=table([seq(op([i,1],prov)=i,i=1..n)]);
    lab2:=table([seq(i=op([i,1],prov),i=1..n)]);
    [seq(table([seq(j=lab[perms[i][lab2[j]]],j=1..n)]),i=1..l)]
end;

groupe:=proc(perms)
# It creates the group
local l,prov,n,lab,lab2,i,j,tours;
    l:=nops(perms);
    prov:=op(op(perms[1]));
    n:=nops(prov);
    lab:=table([seq(op([i,1],prov)=i,i=1..n)]);
    lab2:=table([seq(i=op([i,1],prov),i=1..n)]);
    tours:=[seq(table([seq(j=lab[perms[i][lab2[j]]],j=1..n)]),i=1..l)];
    permgroup(n,{seq(convert([seq(tours[i][j],j=1..n)],'disjcyc'),i=1..l)});
end;