#This maple session computes the functions occuring in the explicit 
#Weierstrass preparation theorem explained in the article
#Rational points on hyperelliptic curves and an explicit 
#Weierstrass preparation theorem in Manuscripta Mathematica.
#By the way, this program bound the number of solutions of systems of 
#two power series in two variables which is required for some cases for
#the elliptic curve Chabauty method.
#
#Main programs are wpt and ellcha2.
#
#For any given power series f in Z_p<n_1,n_2> modulo p^prec, wpt returns
#a polynomial g in n_1 with coefficients in Z_p<n_2> such that f=g*h
#where h is a unit in Z_p<n_1,n_2> (and so f(n_1,n_2)=0 implies that
#g(n_1,n_2)=0). The polynomial g is given modulo p^prec as a vector
#(the i-th component is the term of degree i-1). Moreover wpt returns
#the new precision to consider. In fact, one of the coefficients of f 
#must be invertible modulo p. If this is not the case, f must be divided 
#by p and f/p is then known modulo p^{prec-1}.
#
#The program ellcha2 bound the number of solutions of systems of two
#power series f1 and f2 in two variables n1 and n2 given with the 
#precision p^prec. The user have to give the known solutions in n_2 with 
#their multiplicity in the resultant of the polynomials g1 and g2 
#(obtained by wpt). If  the multiplicity is not known, the program
#will help you. The first result is the number of solution in n2 and 
#the second result is the number of solution in n1 for each known 
#solution in n2 given in the entries.


with(linalg):
#We first want to compute the fonction hi in terms of the functions gi.
#For this we use proposition 1. In this proposition, we need to 
#compute a sum over s-uplet i_0,..,i_{s-1} such that i_0+...+i_{s-1}=k. 
#We represent those s-uplets by vectors.
#compnonnul gives the first non zero component of a vector v
compnonnul:=proc(v)    
local i,s,t;          
s:=vectdim(v):
t:=1:
for i from 1 to s do
    if v[i]=0 then t:=t+1 else break 
fi: od:
if t=s+1 then print("zero vector") fi:
t;
end:

#For any given vector v (with the sum of its coordinates equal to k), 
#suivant gives vectors of the same form for which the first non zero 
#component has been reduced of 1
suivant:=proc(v)
local W,c,i,s,w;
W:={}:
c:=compnonnul(v):
s:=vectdim(v):
for i from c+1 to s do
    w:=v:
    w[c]:=v[c]-1:
    w[i]:=v[i]+1:
    W:=W union {w}:
od:
W:
end:

#finally partition gives all the required s-uplet (the partitions of k)
partition:=proc(s,k)
local vi,vf,i,S ;
vi:=[seq(0,i=1..s)]:
vf:=[seq(0,i=1..s)]:
vi[1]:=k:
vf[s]:=k:
S:={vi}:
while member(vf,S)=false do
      for i from 1 to nops(S) do
          S:=S union suivant(S[i]):
      od:
od:
S;
end:

#indice computes the index of the fi occuring in the formula
indice:=proc(n,k,s,i)
local ind;
ind:=n+s+sum('(s-(j-1))*i[j]','j'=1..s);
ind;
end:

#multinomial computes the multinomial coefficient occuring in 
#(g_0+...+g_{s-1})^k
multinomial:=proc(i,k)
local c,s;
s:=vectdim(i):
c:=k!/product('i[j]!','j'=1..s):
c;
end:

#sommepart computes the finite sum for any given k occuring in the 
#formula
sommepart:=proc(f,g,k,n,p,prec)
local sum,s,S,j,i,a,gi,find; 
sum:=0:
s:=vectdim(g)-1:
S:=partition(s,k):
for j from 1 to nops(S) do
    i:=S[j]:
    find:=coeff(f,n1,indice(n,k,s,i)) mod p^prec:
    gi:=vector(s,0):
    for a from 1 to s do 
        if g[a]=0 and i[a]=0 then 
           gi[a]:=1 else
           gi[a]:=g[a]^(i[a]) mod     p^prec: 
        fi:
    od: 
    sum:=sum+multinomial(i,k)*find*product('gi[l]','l'=1..s) mod p^prec:
od:
sum;
end:
 
#inversemodp computes the inverse of a polynomial in n2 modulo p with
#the precision p^prec
inversemodp:=proc(Pol,p,prec)
local P,Pi,a0,a0i,k;
P:=Pol mod p^prec:
a0:=coeff(P,n2,0):
a0i:=1/a0 mod p^prec:
P:=P*a0i mod p^prec:
if ((P-1) mod p) = 0 then 
                     else print(Pol,"is not invertible modulo", p^prec):
fi:
Pi:=simplify(a0i*(1+sum('(1-P)^k','k'=1..prec-1))) mod p^prec:
Pi;
end:

#Finally, calculhn computes the function hn with the precision p^prec if
#f and g are given with the precision p^{prec-1}
calculhn:=proc(f,g,h,n,p,prec)     
local kmax,s,indmax,hn,k,F;         
s:=vectdim(g)-1:
F:=f mod p^prec:
indmax:=degree(f,n1):
kmax:=max(indmax-s-n,0):
hn:=sum('expand((-1)^k *inversemodp(expand(g[s+1]^(k+1)) mod p^prec,
        p,prec)*sommepart(F,g,k,n,p,prec))' ,'k'=0..kmax) mod p^prec:
hn;
end:

#The proposition 1 is now entirely computable, so that we can now 
#apply Weierstrass prepareation theorem. explwpt return the polynomial
#g and the function h with the precision p^prec if f is general of 
#order s and given with the same precision. explwpt also return test 
#which test if f=g*h modulo p^prec
explwpt:=proc(f,s,p,prec)
local F,i,g,h,k0,j,l,n,indmax,H,G,test;
F:=f mod p^prec:
g:=vector(s+1,0):
indmax:=degree(F,n1):
h:=vector(max(s+1,s+indmax-2*s+1)+1,0):
for i from 1 to s+1 do
    g[i]:=coeff(F,n1,i-1) mod p:
    h[i]:=0:
od:
h[1]:=1:
for k0 from 2 to prec do
    for j from 2 to s+1 do 
        
        h[j]:=sort(calculhn(F,g,h,j-1,p,k0),n2):
        print(h[j]);
    od:
    for j from 1 to s+1 do
    g[j]:=sort(expand(coeff(F,n1,j-1)-sum('h[l]*g[j-l+1]','l'=2..j)), n2)
                                                                mod p^k0:
    od:
od:
for n from s+1 to s+indmax-2*s+1 do
    h[n+1]:=sort(calculhn(F,g,h,n,p,prec),n2):
od:
H:=expand(sum('h[i]*n1^(i-1)','i'=1..vectdim(h))) mod p^prec:
G:=expand(sum('g[i]*n1^(i-1)','i'=1..s+1)) mod p^prec:
test:=expand(F-G*H) mod p^prec:
[test,evalm(g),evalm(h)]:
end:

#general test if f (given with the precision p^prec) is general in n1 
#and give the order and the new precision to consider (in fact, if f 
#is divisible by p, this program compute this division so that we 
#lose in precision)
general:=proc(f,p,prec)
local F,s1,t,s,newprec;
F:=f:
newprec:=prec:
while (F mod p)=0 do
       F:=simplify(F/p):
       newprec:=newprec-1:
od:
s1:=degree(subs(n2=0,F) mod p,n1):
if s1=-infinity  then print(f,"is not general in",n1,"try a linear 
combination with the other power series or exchange n1 and n2 (in the 
two power series)"): break: fi:
[F,newprec,s1]:
end:

#wpt is the weierstrass preparation theorem : given a power series f in
#n1 and n2 with the precision p^prec, wpt apply weierstrass preparation 
#theorem for f relatively to n1. It returns g (h is useless but wpt test
#if f=g*h) as a vector with the new precision to consider (if f is 
#divisible by p, we lose in precision)
wpt:=proc(f,p,prec)
local F,newprec,s,t,test,g,h,A,B;
A:=general(f,p,prec):
F:=A[1]:newprec:=A[2]:s:=A[3]:
B:=explwpt(F,s,p,newprec):
test:=B[1]:g:=B[2]:h:=B[3]:
if test=0 then else
print("bug, please send your session to duquesne@math.u-bordeaux.fr, 
       many thanks."):break: 
fi:
[evalm(g),newprec]:
end:

#Finally, ellcha2 bound the number of solutions of systems of two power 
#series f1 and f2 in two variables n1 and n2 given with the precision 
#p^prec. The user have to give the known solutions of the resultant of
#the polynomial g1 and g2 with their multiplicity. But if this is not 
#the case, the program will help you. The first result is the number 
#of solutions in n2 and the second result is the number of solutions
#in n1 for each known solution in n2 given in the entries.
ellcha2:=proc(f1,f2,p,prec,sol)
local g1,g2,G1,G2,G1u,G2u,res,newprec,s,S,i,F,sol2,s2,lowdeg;
g1:=wpt(f1,p,prec):
g2:=wpt(f2,p,prec):
newprec:=min(g1[2],g2[2]):
print('modulo',p^newprec):
G1:=sort(sum(g1[1][i]*n1^(i-1),i=1..vectdim(g1[1])) mod p^newprec,n1): 
print("first polynom",G1):
G2:=sort(sum(g2[1][i]*n1^(i-1),i=1..vectdim(g2[1])) mod p^newprec,n1):
print("second polynom",G2):
res:=sort(resultant(G1,G2,n1) mod p^newprec,n2):
print("resultant",res):
if (res mod p^newprec)=0 then print("precision too low or resultant 
                                     equals to zero") else
  while (res mod p)=0 do
       res:=simplify(res/p):
  od:
  s:=degree(res mod p,n2):
  print("there is at most",s,"solutions in n2"):
  if (vectdim(sol)=s and s<>0) then
     sol2:={}:

     for i from 1 to s do
         if member(sol[i],sol2)=false then sol2:=(sol2 union {sol[i]}):
         fi:
     od:
 s2:=nops(sol2):
     S:=vector(s2,0):
     for i from 1 to s2 do
         F:=subs(n2=sol2[i],f1):
         if (F mod p^newprec)=0 then print("precision too low") else
            while (F mod p)=0 do
              F:=simplify(F/p):
            od:
            S[i]:=degree(F mod p,n1):
            print("for n2 =",sol[i],"there is at most",S[i],"solutions 
                  in n1 corresponding to the zeros of",F):
         fi:
     od:
  else 
    lowdeg:=ldegree(res):
    if ldegree=0 then     
       print("You have to think to indicate all the known solutions for
             n2. If done, this means that elliptic curve Chabauty does 
             not work in this case"):break: 
    else print("Be careful, 0 can be a multiple zero. In this case it is
               probably of order",lowdeg,"for the resultant. You must in
               dicate it in the entries of the program. For exemple, for
               the point at infinity, 0 is always of order 4"):break:
    fi:
  fi:
fi:
[s,evalm(S)]:
end: