read "puiseuxmodp.mpl": read "choixp.mpl": mypuiseux:=proc(f,x,y,pol,N) # This algorithm computes numerical evaluation of Puiseux expansions # above the roots of pol. # f is a polynomial in x and y. The polynomial is supposed to be known # exactly. # pol is a polynomial in x. It is supposed sqrfree. NB: if we know # that the polynomail is irreducible, than we only need to compute the # polygone tree above one root of pol modulo p. Thus, this case would # need a specific treatment that is not done here (we suppose that # such information has a too big cost, and that we avoid # factorisations, and so irreducibility criterion, in a field of # characteristic 0). # N is the number of terms that we want to compute. N=0 if we want the # singular part. # Eight options are available: # T an unknown which is the parameter (in this case, the algorithm # output the Puiseux expansions expressed in parametrization) # valuation=c gives the valuation c in pol 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. # primenumber=p gives the primenumber that we want to use during the # modular algorithm # probability=eps gives the probability we use for the Monte Carlo # Good Prime algorithm # trailingcoefficient=tcRF gives the trailing coefficient of # resultant_y(F,F_y)for the Las Vegas Good Prime algorithm # roots = numroots gives the numerical roots of pol. # 'moreterms' if we want to know the data c such that # G(x,y)=F(P,Q)/x^c to compute more terms later (but we just compute # the singular part here) # 'irreducible' if the polynomial p is irreducible over the field we # are considering (this is to improve the modular part of the algorithm) # Remark: one input should be the precision that we require for the # coefficients of the power series, or something like that (for # instance, the precision of the input). This is not the cas at this time. local opt,c,p,eps,tcRF,t,T,i,trees,F,res,ext,K,M,prov,dy,dx,numF,co,numroot,numpols,z,other,j,more,irr,t0,t1,t2; # We first consider the option given in input. t1:=time(); userinfo(1,'puiseux',`Computing Puiseux expansions above roots of`,pol); `if`(N=0,userinfo(1,'puiseux',`We are searching the singular parts`),userinfo(1,'puiseux',`we want`,N,`terms`)); opt:=[args[6..nargs]]; for i in opt do if op(1,i)='valuation' then c:=op(2,i) elif op(1,i)='primenumber' then p:=op(2,i) elif op(1,i)='probability' then eps:=op(2,i) elif op(1,i)='trailingcoefficient' then tcRF:=op(2,i) elif op(1,i)='roots' then numroot:=op(2,i) elif i='moreterms' then more:=true elif i='irreducible' then irr:=true elif i::name then t:=i else error "wrong input" fi; od; if not assigned(more) then more:=false fi; if not assigned(irr) then irr:=false fi; # We want only one primitive element for the coefficient field. Thus, # we begin to create this primitive element if necessary. t0:=time(); userinfo(4,'puiseux',`working on extensions`); 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; F:=f; fi; userinfo(4,'puiseux',`running time`,time()-t0); t0:=time(); userinfo(4,'puiseux',`finding a good prime number`); if not assigned(p) then # For the moment, we suppose that if p is not given in the input, then # pol = x. Otherwise, we need to perform the change of variable # x <- x + RootOf(pol) on F to get a local criteria. # For a final version, the hypothesis should only be that pol is # irreducible... if assigned(eps) then p:=MCGP(F,x,y,M,z,eps) elif assigned(tcRF) then p:=LVGP(F,x,y,M,z,tcRF) else # In this last case, I decide to use a Las Vegas strategy, and so to # compute the resultant... tcRF:=tcoeff(evala(eval(resultant(F,diff(F,y),y),x=x+RootOf(pol)))); p:=LVGP(F,x,y,M,z,tcRF) fi; fi; # p is the prime number that we will use in the modular algorithm. userinfo(4,'puiseux',`running time`,time()-t0); userinfo(2,'puiseux',`prime number used will be`,p); if not assigned(t) then t:=T fi; t0:=time(); userinfo(2,'puiseux',`Computing the polygon trees`); c,trees:=computetrees(F,x,y,M,z,p,pol,irr,`if`(assigned(c),c,NULL)); userinfo(2,'puiseux',`running time`,time()-t0); t2:=time(); userinfo(2,'puiseux',`Beginning numerical computations`); dy:=degree(f,y); dx:=degree(f,x); numF:=expand(evalf(F)); co:=op(map(abs,[coeffs(numF,{x,y})])); #ceil(Digits+log10(max(co))-log10(min(co))+log10(dy^3*dx^2)); Digits:=ceil(Digits+evalf(log10(max(co))-log10(min(co))+log10(dy^3*dx^2))); # for the precision : we take care of the size of the coefficients of # F and we add some Digits according to the complecity result. # We recall that this implementaion is not a certified one userinfo(2,'puiseux',`Digits used`,Digits); numF:=expand(evalf(F)); if not assigned(numroot) then t0:=time(); userinfo(3,'puiseux',`Computing roots of the discriminant factor`); numroot:=[fsolve(pol,x,'complex')]; userinfo(3,'puiseux',`running time`,time()-t0); fi; t0:=time(); userinfo(3,'puiseux',`Computing numerical shifts`); numpols:=[seq([polytoArray(numshift(numF,x,i),x,y,dy,min(dx,c)),t,y,0,i],i in numroot)]; userinfo(3,'puiseux',`running time`,time()-t0); t0:=time(); userinfo(3,'puiseux',`Computing singular parts`); prov:=mynewton(numpols,trees,N+c,t,y,evalb(N=0)); userinfo(3,'puiseux',`running time`,time()-t0); res:=NULL; if N=0 then prov:=[eval(prov,y=0)]; for i in numroot do res:=res,[seq(`if`(j[4]=i,[x=j[1]+i,y=j[2],`if`(more,j[3],NULL)],NULL),j=prov)]; od; res:=[res]; else t0:=time(); userinfo(3,'puiseux',`Computing more terms`); for i in numroot do res:=res,[seq(`if`(j[5]=i,[Arraytopoly(j[1],x,y),j[2]+i,j[3]],NULL),j=prov)]; od; res:=[res]; for i to nops(res) do for j to nops(res[i]) do other:=regular(res[i,j,1],diff(res[i,j,1],y),x,y,0,N+1-degree(lcoeff(res[i,j,3],y),t),true); res[i,j]:=[x=res[i,j,2],y=eval(eval(res[i,j,3],y=convert(other,'horner')),x=t)]; od; od; userinfo(3,'puiseux',`running time`,time()-t0); fi; userinfo(2,'puiseux',`numerical computations running time`,time()-t2); userinfo(1,'puiseux',`total running time`,time()-t1); res; end; computetrees:=proc(F,x,y,M,z,p,pol,irr) # irr is true if pol is known irreducible and false otherwise. local xp,extp,rootp,trees,i,Fp,prov,c,degfac,k,mini; if nargs=9 then c:=args[9] fi; xp,extp:=oneroot(M,z,p,[]); # xp is the root of pol in Fp which generate the smallest extension. # now we have to consider the factorisation of pol in Fp(extp) rootp:=myRoots(eval(pol,RootOf(M)=xp),x,p,extp); # rootp 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 # D is op(extp) # E is the expression of D in Fp(C) if irr then degfac:=seq(degree(eval(op(i[3]),_Z=z),z),i in rootp); mini:=min(degfac); member(mini,[degfac],'k'); rootp:=[rootp[k]]; userinfo(2,'puiseux',`We will work in an extension of degree`,mini); fi; trees:=NULL; for i in rootp do if i[2]<>1 then error "bad prime p" fi; # In such a case, the square free factorisation has changed. Thus, the # polygon trees over some critical points have changed. Fp:=Expand(eval(eval(F,RootOf(M)=xp),i[4]=i[5])) mod p; # We write the reduction of F in the new extension prov:=puiseuxmodp(Fp,x,y,i[1],0,p,`if`(i[3]=[],[],[i[3]]),`if`(assigned(c),c,NULL)); trees:=trees,op(1,prov)$`if`(irr,degree(pol,x),degree(eval(op(i[3]),_Z=z),z)/degree(eval(op(i[4]),_Z=z),z)); if not assigned(c) then c:=op(2,prov) fi; od; c,[trees] end; oneroot:=proc(poly,x,p,alpha) # It computes one of the roots of poly (poly is supposed to belong to # F_p(alpha)[x]) in the algebraic closure of Fp which adds the # smallest extension. The output gives also the extension of F_p # where this root belongs. # Here, alpha is [] or [RootOf()]. local beta,fac,degfac,mini,i,k,sol,delta; beta:=op(alpha); fac:=Factors(poly,beta) mod p; degfac:=seq(degree(fac[2,i,1],x),i=1..nops(fac[2])); mini:=min(degfac); member(mini,[degfac],'k'); sol:=fac[2,k,1]; if mini > 1 then if beta=NULL then beta:=RootOf(sol); sol:=beta; else delta:=newext(beta,RootOf(sol),p); sol:=Roots(eval(sol,beta=Roots(eval(op(beta),_Z=x),delta) mod p),delta) mod p # Here, we express RootOf(sol) in Fp(delta) fi; else sol:=-coeff(sol,x,0) mod p fi; sol,[beta] end; numshift:=proc(P,x,x0) expand(PolynomialTools['Translate'](evalf(P),x,x0)) end; evaldiagnum:=proc(F,Nc,m,q,xi,l) # This function computes F(x^q,x^m(xi^(1/q)+y))/x^l where F is # represented by a table. # xi and the elements F[i,j] are supposed to be floating points # numbers. local G,k,d,P,i,y,dX,dY,N,a,t; dY:=op([1,2],[ArrayDims(F)]); 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 slopes do not necesseraly increase). 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 and (k-m*i)/q >=0 ,F[i,(k-m*i)/q]*y^i,0 ),i=0..d); P:=numshift(P,y,evalf(xi)^(1/q),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: mynewton:=proc(numpols,trees,N,T,y,sing) # This recursive algorithm computes numerical approximations of # singular parts of Puiseux series of polynomials in numpols above # 0. If our im is to compute n terms or less of each series, then N # equals n plus a bound. This is for the truncations of the power of X # during the algorithm (see lemma 16 of my thesis) # trees is a list of polygon trees # numpols is a list of [H,X,Y] where H is a numerical polynomial # expressed in a table, and X and Y give the parametrization # associated to the polynomial H. # dy and c give the size of the tables. # sing is a boolean which is true if we are just interested by the # singulart part of the Puiseux expansions and false otherwise. # t and y are the parameter in which we want to express the # parametrizations. local res,sameN,Np,i,samepart,roots,j,k,Delta,m,q,l,mult,alpha,beta,H,X,Y,delta,xi,Ri,Li,s,n,h,x0,z,c; res:=NULL; sameN:=[tria(numpols,trees,N)]; for i in sameN do Np:=i[2,1]['polygon']; # Np is the common Newton polygon of polynomials in i[2]. samepart,roots:=trib(Np,op(i)); for k to nops(Np)-1 do Delta:=Np[k..k+1]; m,q,l:=op([1,1..3],i[2,1]['slope',k]); for j in samepart do # j = [P,L] where P is a list of [H,X,Y] and L is the list of polygon # trees associated. All these tree have the same Newton polygon Np as # root vertex, and each characteristic polynomial associated to each # edge of Np has the same multiplicity structure. mult:={seq(h[2],h=j[2,1]['slope',k][2])}; # mult is the set of multiplicities of the root of the characteristic # polynomials associated to Delta. for alpha in mult do if alpha<>1 then Ri:=NULL; # R1 will store all the ne polynomials defined by the edge Delta and # each root of multiplicity alpha of the phi_Delta. for beta in j[1] do H,X,Y,c,x0:=op(beta); #print("coucou",Arraytopoly(H,'x','y')); xi:=roots[H,Delta]; Ri:=Ri,seq([evaldiagnum(H,N,m,q,delta,l),eval(X,T=T^q),eval(Y,{T=T^q,y=T^m*(y+delta^(1/q))}),q*c+l,x0],delta in xi[alpha,1]); od; Ri:=[Ri]; Li:=NULL; # Li will store all the associated subtrees. for beta in j[2] do s:=beta['slope',k][2]; Li:=Li,seq(`if`(s[n,2]=alpha,beta['rec',k,n]$(degree(eval(op(s[n,3]),_Z=z),z)/degree(eval(op(s[n,4]),_Z=z),z)),NULL),n=1..nops(s)); od; Li:=[Li]; res:=res,mynewton(Ri,Li,N,T,y,sing) else for beta in j[1] do H,X,Y,c,x0:=op(beta); xi:=roots[H,Delta]; res:=res,seq([`if`(sing,NULL,evaldiagnum(H,N,m,q,delta,l)),eval(X,T=T^q),eval(Y,{T=T^q,y=T^m*(y+delta^(1/q))}),q*c+l,x0],delta in xi[alpha,1]); od; fi; od; od; od; od; res; end; tria:=proc(numpols,trees,c) # This function sort the polynoms according to their generic (or # excptionnal) Newton polygons. local polygons,i,M,j,belongs,notbelongs,pt,x,y,maxval,val,zero,mini,minpol,other; polygons:=seq(i['polygon'],i=trees); if nops({polygons})=1 then return([numpols,trees]) fi; polygons:=[polygons]; M:=op([1,nops(polygons[1]),1],polygons); # M is the width of the polygons minpol:=polygon(polytoArray(add(add(x^j[2]*y^j[1],j in i),i in polygons),x,y,M,c),M); # minpol is the lower part of the convex hull of the union of the # Newton polygons in polygons. for i in minpol do belongs:=NULL; notbelongs:=NULL; for j to nops(polygons) do if i in polygons[j] then belongs:=belongs,j; else notbelongs:=notbelongs,j; fi; od; notbelongs:=[notbelongs]; if notbelongs <> [] then pt:=i; break; fi; od; # Here we have one point which belongs to some Newton polygons but is # "under" the others. val:=[seq(abs(i[1][pt[1],pt[2]]),i in numpols)]; maxval:=max(op(val)); zero:=NULL; for i to nops(notbelongs) do mini:=min(op(val)); member(mini,val,'j'); val[j]:=maxval+1; zero:=zero,j; od; zero:={zero}; other:={$1..nops(val)} minus zero; zero:=[seq(numpols[i],i in zero)]; other:=[seq(numpols[i],i in other)]; belongs:=[seq(trees[j],j in [belongs])]; notbelongs:=[seq(trees[j],j in notbelongs)]; tria(zero,notbelongs,c),tria(other,belongs,c) end; trib:=proc(N,numpol,trees) # It separates polynomials according to their characteristic # polynomials and the associate partitions. local R,i,Delta,d,R1,j,polcar,k,l,partitions,R0,Z,S,v,xi,m; S:=table(); R:=[[numpol,trees]]; for i to nops(N)-1 do Delta:=N[i..i+1]; d:=(Delta[2,1]-Delta[1,1])/trees[1]['slope',i][1,2]; R1:=NULL; for j in R do polcar:=[seq([[[phi(k[1],op(Delta),Z),1,1]],k],k in j[1])]; partitions:=[seq([[seq(l[2]$(`if`(l[3]=[],1,degree(eval(op(l[3]),_Z=Z),Z))/`if`(l[4]=[],1,degree(eval(op(l[4]),_Z=Z),Z))),l=k['slope',i][2])],k],k in j[2])]; R0:=[tric(polcar,partitions,d,Z)]; R1:=R1,seq([[seq(k[2],k in l[1])],[seq(k[2],k in l[2])]],l in R0); for m in R0 do for k in m[1] do # We now compute roots of the characteristic polynomial using the # list of gcd triplets. v:=[seq(l[2], l in k[1])]; xi:=[seq([[fsolve(myagcd(v[-l-1],v[-l-2],Z,degree(v[-l-2],Z))[2],Z,'complex')],l],l=1..nops(v)-2),[[fsolve(v[1],Z,'complex')],nops(v)-1]]; S[k[2,1],Delta]:=xi; od; od; od; R:=[R1]; od; R,S; end; tric:=proc(triplets,partitions,d,Z) # It separates polynomials according to the gcd's of (phi,phi') and # other successives gcd's # triplets is the list of [U,H] where U is a list which gives the r # first gcd triplets associated to the characteristic polynomial # defined by the bivariate polynomial H. These triplets are ordered # such that the last gcd computed is the first of the list. # partitions is the list of partitions associated to the # characteristic polynomials (each one is given with its associated # tree) # d is the degree of each characteristic polynomial # Z is the variable in which gcd triplets are expressed local n,degs,i,m,singval,mins,N,ou1,ou2,j,R,L1,L2,prov; if d=0 then return([triplets,partitions]) fi; R:=NULL; n:=nops(triplets); degs:=[seq(d-nops(i[1]),i in partitions)]; # degs is the list of the degrees for gcd(phi,phi') m:=add(i,i in degs); # m is the number of singular values which are symbolically equal to # zero singval:=[seq(sort(evalf(LinearAlgebra['SingularValues'](LinearAlgebra['SylvesterMatrix'](i[1,1,1],diff(i[1,1,1],Z),Z),'output'='list')),`<`),i in triplets)]; mins:=[seq(i[1],i in singval)]; N:=[0$n]; while m>0 do member(min(op(mins)),mins,'i'); N[i]:=N[i]+1; mins[i]:=`if`(N[i]=d,infinity,singval[i,N[i]+1]); m:=m-1; od; for i in {op(N)} do ou1:=NULL; while member(i,N,'j') do ou1:=ou1,j; N[j]:=-1; od; ou1:=[ou1]; # ou1 contains the indices of the characteristic polynomials for which # degree(gcd(phi,phi'))=i ou2:=NULL; while member(i,degs,'j') do ou2:=ou2,j; degs[j]:=-1; od; ou2:=[ou2]; # ou2 contains the indices of trees which have a partition l such that # d-cardinal(l)=i if nops(ou1) <> nops(ou2) then error "symbolic and numeric values do not correspond" fi; # At this stage, we have separate the polynomials according to the # degrees of gcd(P,P'). To sort them according to the multiplicities # structure, we need recursive calls. But we first have to compute the # approximate gcd. L1:=NULL; for j in ou1 do prov:=triplets[j]; prov[1]:=[myagcd(prov[1,1,1],diff(prov[1,1,1],Z),Z,i),op(prov[1])]; L1:=L1,prov; od; L1:=[L1]; L2:=[seq([diminuer(partitions[j,1]),partitions[j,2]],j in ou2)]; R:=R,tric(L1,L2,i,Z); od; R; end; diminuer:=proc(l) local i; [seq(`if`(i>1,i-1,NULL),i in l)] end; myagcd:=proc(p,q,x,m) # It computes the approximate gcd triplet of P and Q in the variable # x, knowing its degree m. local P,Q,cP,cQ,dp,dq,S,V,u,v,w,i,l; if m=0 then return [1,p,q] fi; P,Q:=evalf(expand(p)),evalf(expand(q)); dp,dq:=degree(P,x),degree(Q,x); cP:=op(convert(PolynomialTools['CoefficientVector'](P,x),list)); cQ:=op(convert(PolynomialTools['CoefficientVector'](Q,x),list)); S:=LinearAlgebra['Transpose'](Matrix([seq([0$i,cP,0$(dq-m-i)],i=0..dq-m),seq([0$i,cQ,0$(dp-m-i)],i=0..dp-m)])); V:=LinearAlgebra['SingularValues'](S,'output'='Vt'); V:=V[dp+dq+2*(1-m),1..(dp+dq+2*(1-m))]; v:=add(V[l+dq-m+1]*x^(l-1),l=1..dp-m+1); # v*u=P where u = gcd(P,P') # maybe a bug with Singular Values (or I do not understand) # The following line is to correct this v:=eval(v,I=-I); w:=add(V[l]*x^(l-1),l=1..dq-m+1); # w*u=P' # maybe a bug with Singular Values (or I do not understand) # The following line is to correct this w:=eval(w,I=-I); u:=quo(P,v,x); # u is the gcd. A better way to compute u from v should be # implemented. This way may not be numerically stable. One solution # should use the "least square division" (see z. zeng papers for # instance). I will consider this point later. [u,v,w]; end; phi:=proc(F,a,b,T) # It computes phi for one vertex Delta=[a,b] local ptI,ptJ,i,j,m,q,l,phi,k; ptI,ptJ,i,j:=`if`(a[1]