polytoArray:=proc(Pf,x,y,dY,N)
# Starting from a polynomial Pf in K[x,y], it creates an Array F such
# that F[i,j] is the coefficient of Pf in x^j*y^i.
local F,a,t,i,dx,dy;
    F:=Array(0..dY,0..min(N,degree(Pf,x)));
    a:=[coeffs(Pf,[x,y],'t')];
    t:=[t];
    for i to nops(t) do
        dy,dx:=degree(t[i],y),degree(t[i],x);
        if dx<=N then F[dy,dx]:=a[i] fi;
    od;
    F;
end;


newRootOf:=proc(beta,alpha,p)
# Given 2 RootOf beta (which does not depend on any other RootOf) and
# alpha (which can depend of beta), it computes the RootOf delta such
# that Fp(delta)=Fp(beta,alpha) and it gives beta and alpha expressed
# in Fp(delta).

# The polynomial which defines beta is supposed to be irreducible in
# Fp, and the one which defines alpha is supposed to be irreducible in
# Fp(beta).
local P,Q,T,Z,z,delta,beta1,alpha1;
option remember;
    P:=eval(op(beta),_Z=T);
    Q:=eval(op(alpha),{beta=T,_Z=Z});
    delta:=RootOf(Nextprime(z^(degree(P,T)*degree(Q,Z)),z) mod p);
# delta=RootOf(R) where R is an irreducible polynomial over Fp. Its
# degree is the degre of the extension Fp(alpha,beta). Thus,
# Fp(alpha,beta)=Fp(delta).
# We now have to compute the polynomials P1 and Q1 such that
# alpha=P1(delta) and beta=Q1(delta)
    beta1:=modp(Roots(P,delta),p)[1,1];
    alpha1:=modp(Roots(eval(Q,T=beta1),delta),p)[1,1];
    [delta,beta1,alpha1]
end;



myRoots:=proc(pol,z,p,alpha)
# This algorithm computes the roots of the polynomial pol(z) in
# F_p(alpha). More precisely, the algorithm make the factorization of
# pol(z) in Fp(alpha)[z] and then output one root of each irreducible
# factor (expressed in the RootOf form). It also gives the
# multiplicity of this root and the extension of Fp in which this root
# belongs (and the root is expressed in this extension).

# The output is a list of list [A,B,C,D,E] where

# A is a solution of pol expressed in Fp(C)

# B is the multiplicity of the root A in pol

# C is a primitive element above Fp or []

# D is op(alpha) or []

# E is the expression of D in Fp(C)

# alpha is [RootOf(...)] or []
local fac,i,sol,r,beta;
    beta:=op(alpha);
    fac:=op(2,Factors(pol,beta) mod p);
    sol:=NULL;
    for i in fac do
        if degree(i[1],z)>1 and type(beta,RootOf) then
            r:=newRootOf(beta,RootOf(i[1]),p);
            sol:=sol,[r[3],i[2],r[1],beta,r[2]];
        elif beta=NULL and degree(i[1],z)>1 then
            sol:=sol,[RootOf(i[1]),i[2],RootOf(i[1]),[],[]];
        elif type(beta,RootOf) then
            sol:=sol,[-coeff(i[1],z,0) mod p,i[2],beta,beta,beta]
        else
            sol:=sol,[-coeff(i[1],z,0) mod p,i[2],[],[],[]]
        fi;
    od;
    [sol];
end;


shift:=proc(P,x,c,p)
# It evaluates the polynomial P viewed as a univariate polynomial in x
# in x+c (in FP(alpha) )
#    Expand(eval(P,x=x+c)) mod p;
# I try fast algorithms :
Expand(PolynomialTools[Translate](P,x,c)) mod p
end;




puiseuxmodp:=proc(F,x,y,pt,N,p,alpha)
# The modular rational Newton-Puiseux algorithm
# It computes the N first terms of Puiseux series of F above pt in the
# algebraic closure of Fp[[x]]. If N=0, then it computes the singular
# part of the Puiseux series.

# F is supposed to be in the two variables x and y and with
# coefficients in F_p(alpha) (here, alpha is [] or [RootOf(...)]).

# pt is supposed to be given in Fp(alpha). Otherwise, some problems
# will appear since Maple does not support more than one algebraic
# extension above Fp.

# Two options are available: T an unknown which is the parameter (in
# this case, the algorithm output the Puiseux expansions expressed in
# parametrization) and c the valuation in (x-pt) of the discriminant
# of F in Y (this data is used to have a faster algorithm but can take
# time to compute). c must be an integer.

# Remark: the regular part is not yet implemented, and so does not
# appear in this version of the code.
local dy,r,res,G,prov,c,T,parametrization,arg;
    for arg in [args[8..nargs]] do
        if arg::integer then
            c:=arg;
        elif arg::name then
            T:=arg;
            parametrization:=true;
        else
            error "bad arguments"
        fi;
    od;
    if eval(F,x=0)=0 then error "x divides F" fi;
    dy:=degree(F,y);
    G:=shift(F,x,pt,p);
    if not assigned(c) then
        c:=ldegree(Expand(resultant(G,diff(G,y),y)) mod p)
    fi;
    if c=0 then return("regular polynomial") fi;
# c is the valuation of the discriminant of G in x. This is an upper
# bound for the number of terms that we will have to compute to get
# the singular part. Thus, we do not have to consider coefficients
# with a power greater than c in x to compte the singular part.
    G:=polytoArray(G,x,y,dy,c);
    prov:=newtonmodp(G,dy,dy,max(c,N),p,alpha,evalb(N=0),true);
# newtonmodp computes the singular part of the series.
    if N=0 then
        res:=prov;
# I have to add the regular part to get the N terms.
    fi;
# Don't forget to add the output in a series form and not in a
# parametrization's one (according the value of parametrization).
    [res,c]
end;



newtonmodp:=proc(F,dY,M,N,p,ext,sing,first)
# Remark: the case where Y divides F is not yet implemented.
# -> done with the generic newton polygon

# This recursive algorithm computes Puiseux series of F above (0,0).

# F is supposed to be in F_p(op(extension))[x,y] but is given in a
# matrix form (this is faster to get access to coefficients).

# dY and N give the size of the table.

# M gives the restriction of terms that we have to consider to
# compute the Newton polygon of F (it is the multiplicity of the
# previous step.

# sing is a boolean which is true if we are just interested by the
# singulart part of the Puiseux expansions and false otherwise.

# first is a boolean which is true if we are at the initial call and
# false if we are at a recursive call

# The output is given in the following way:
# res['polygon'] will give the Newton polygon of F.
# res['slope',i] gives the number m,q,l,u,v associated to the i-th edge of
# polygon, and the roots of its characteristic polynomial.
# res['rec',i,j] gives the result for the recursive call of the algorithm
# for the i-th edge and the j-th root of the characteristic polynomial.
local res,a,i,Delta,m,q,l,phi,Z,xi,G,j,u,v,newext;
    a:=polygon(F,M,`if`(first,'exceptionnel','generic'));
    res['polygon']:=a;
    for i to nops(a)-1 do
        Delta:=a[i..i+1];
# Delta is the edge of the newtonpolygon we are looking for.
        m,q,l,u,v,phi:=pts(F,op(Delta),Z);
# m,q and l represent the edge Delta.
# phi is the characteristic polynomial asociated to Delta.
        xi:=myRoots(phi,Z,p,ext);
# xi is the list of roots of phi given with multiplicities. More
# precisely, xi[i] is a list [root,multiplicity,delta,alpha,alpha1]
# where delta is the new extension (delta is a RootOf of an
# irreducible polynomial of Fp or []), alpha is the previous extension
# and  alpha1 is alpha1 expressed in a polynomial of delta (alpha=[]
# if there was not any previous extension).
        res['slope',i]:=[[m,q,l,u,v],xi];
# We stock informations which gives the relation between each
# RootOf. It can be useful if we want to express the result with a
# parametrization (as in puiseuxparametrization).
        for j to nops(xi) do
            newext:=`if`(xi[j,3]=[],[],[xi[j,3]]);
            if xi[j,2]<>1 then
                G:=evaldiagmodp(F,dY,N,m,q,u,v,xi[j,1],l,p,ext,xi[j,4]);
                res['rec',i,j]:=newtonmodp(G,dY,xi[j,2],N-m,p,newext,sing,false);
# We can reduce the number of terms that we search by m since this version
# of the algorithm makes factorization on the result (in fact, we
# could reduce by e*m/q but we do not know e)
            elif not(sing) then
                res['rec',i,j]:=evaldiagmodp(F,dY,N,m,q,u,v,xi[j,1],l,p,ext,xi[j,4]);
# Here, we do not compute the new polynomial if we just want the
# singular part of the Puiseux expansion. Otherwise, we have to
# compute it to apply the regular algorithm to the series.
            fi;
        od;
    od;
    res
end;


polygon:=proc(F,M)
# It computes the generic newton polygon associated to the polynom F expressed
# in a table.

# M is the multiplicity of the previous root (or the degree of F in Y
# at the beginning).

# one option: args[3] belongs to {exceptionnel,generic} if we want to
# compute one of this Newton polygon (and not the classic one)
local ii,i,j,a,good,k,N,notseen,pol,which;
    if nargs=3 then which:=args[3] fi;
    N:=op([2,2],[ArrayDims(F)]);
    ii:=Array(0..M);
    for i from 0 to M do
        notseen:=true;
        for j from 0 to `if`(which='generic',min(N,M-i),N) do
# can be min(N,M-i) to reduce the number of coefficients checked
# when we are considering the generic polygon.
            if F[i,j]<>0 then ii[i]:=j; notseen:=false; break fi;
        od;
        if notseen then ii[i]:=N+1 fi;
    od;
    if which='exceptionnel' then ii[0]:=0 fi;
# see the definition of the exceptional polygon
    a:=[seq([i,ii[i]],i=0..M)];
    good:=Array([seq(true,i=1..M+1)]);
    for i from 2 to M do
        for j to i-1 do
            for k from i+1 to nops(a) do
               if a[i][2]>=(a[j][2]-a[k][2])/(a[j][1]-a[k][1])*(a[i][1]-a[k][1])+a[k][2] then
                    good[i]:=false;
                    break;
                fi;
            od;
            if not(good[i]) then break; fi;
        od;
    od;
    pol:=[seq(`if`(good[i],a[i],NULL),i=1..M+1)];
    if which='generic' then
# here we have the classical Newton polygon. We now compute the
# generic Newton polygon.
        for i from nops(pol) to 2 by -1 do
            if pol[i-1,2]>=pol[i,2]+pol[i,1]-pol[i-1,1] then
                pol:=[[0,pol[i,2]+pol[i,1]],op(pol[i..-1])];
                break
            fi;
        od;
    fi;
    pol
end:


pts:=proc(F,a,b,T)
# It computes m,q,l,phi for one vertex Delta=[a,b]
local ptI,ptJ,i,j,m,q,l,phi,k,u,v;
    ptI,ptJ,i,j:=`if`(a[1]<b[1],op([op(a),op(b)]),op([op(b),op(a)]));
    if ptJ=j then
        m,q:=0,1;
    else
        m,q:=-numer((ptJ-j)/(ptI-i)),denom((ptJ-j)/(ptI-i));
    fi;
    if q=1 then
        u,v:=1,0
    else
        igcdex(q,m,'u','v');
        u,v:=u,-v;
    fi;
    l:=m*ptI+q*ptJ;
    phi:=0;
    for k from ptI to i by q do
        phi:=phi+F[k,(l-m*k)/q]*T^((k-ptI)/q)
    od;
    m,q,l,u,v,phi;
end:


evaldiagmodp:=proc(F,dY,Nc,m,q,u,v,xi,l,p,ext,exttonewext)
# This function computes F(xi^v*x^q,xi^u*x^m(1+y)) mod p where F is
# represented by a table.
# The elements F[i,j] are supposed to be in Fp(ext). xi is an element
# of Fp(newext) and exttonewext is ext expressed in a polynoiaml of
# newext. The elements G[i,j] of the output will be expressed in
# Fp(newext).
local G,k,d,P,i,y,dX,N,a,t;
    dX:=op([2,2],[ArrayDims(F)]);
    N:=min(Nc,q*(q*dX+m*dY)-l);
# N is the minimum between the degree in X of the polynomial we are
# computing and the bound that we have for the number of terms.
    G:=Array(0..dY,0..N);
# Here, we cannot reduce the width of the array less than dY since we
# do not know that will be the next slope (we use Duval's algorithm,
# so the slope does not necesseraly increases).
    for k from l to N+l do
# At the next step, we need the monomials up to x^N. Thus, we do not
# have to compute the shift for k=l..q*N, since the values between N+l
# and q*N will generate monomials with degree in x strictly greater
# than N. We have to go up to N+l since we will divide the polynomial
# by x^l.
        d:=`if`(m<>0,min(floor(k/m),dY),dY);
        P:=add(`if`(((k-m*i)/q)::integer and (k-m*i)/q<=dX ,F[i,(k-m*i)/q]*xi^((k-m*i)*v/q)*y^i , 0 ),i=0..d);
        if ext<>[] then
            P:=Expand(eval(P,op(ext)=exttonewext)) mod p
        fi;
        P:=shift(P,y,xi^u,p);
        a:=[coeffs(P,y,'t')];
        t:=[t];
        for i to nops(t) do
            G[degree(t[i],y),k-l]:=a[i]
        od;
    od;
    G
end:


puiseuxparametrization:=proc(L,x,y,T,Y,p)
# This algorithm compute a parametrization in Fp(alpha) of a Puiseux
# series given in a tree form L.

# L is the tree of the Newton polygons and characteristic
# polynomials.

# x and y are the two functions of the parametrization (the first call
# of this algorithm is with x=T and y=Y)

# T is the parameter and Y is the second variable (which is put to 0
# at the end of the algorithm).

# p is a prime number.
local allser,m,q,l,u,v,xi,i,j;
    allser:=NULL;
    for i to nops(L['polygon'])-1 do
        m,q,l,u,v:=op(L['slope',i][1]);
        xi:=L['slope',i][2];
        for j to nops(xi) do
            if xi[j,2]=1 then
                if xi[j,3]<>xi[j,4] and xi[j,4]<>[] then
                    allser:=allser,[Expand(eval(eval(x,xi[j,4]=xi[j,5]),T=xi[j,1]^v*T^q),xi[j,3]) mod p,Expand(eval(eval(y,xi[j,4]=xi[j,5]),{T=xi[j,1]^v*T^q,Y=xi[j,1]^u*T^m}),xi[j,3]) mod p]
                else
                    allser:=allser,[Expand(eval(x,T=xi[j,1]^v*T^q)) mod p,Expand(eval(y,{T=xi[j,1]^v*T^q,Y=xi[j,1]^u*T^m})) mod p]
                fi;
            else
                if xi[j,3]<>xi[j,4] and xi[j,4]<>[] then
                    allser:=allser,puiseuxparametrization(L['rec',i,j],Expand(eval(eval(x,xi[j,4]=xi[j,5]),T=xi[j,1]^v*T^q),xi[j,3]) mod p,Expand(eval(eval(y,xi[j,4]=xi[j,5]),{T=xi[j,1]^v*T^q,Y=xi[j,1]^u*T^m*(1+Y)}),xi[j,3]) mod p,T,Y,p)
                else
                    allser:=allser,puiseuxparametrization(L['rec',i,j],Expand(eval(x,T=xi[j,1]^v*T^q)) mod p,Expand(eval(y,{T=xi[j,1]^v*T^q,Y=xi[j,1]^u*T^m*(1+Y)})) mod p,T,Y,p)
                fi;
            fi;
        od;
    od;
    allser;
end;