# $Source: /u/maple/research/lib/algcurves/src/RCS/puiseux,v $
# $Notify: mvanhoei@daisy.uwaterloo.ca $

macro(
	RAD={radical,nonreal},
	puiseux=`algcurves/puiseux`,
	v_ext_m=`algcurves/v_ext_m`,
	lift_exp=`algcurves/lift_exp`,
	lift_exp_m1=`algcurves/lift_exp_m1`,
	truncate_subs=`algcurves/truncate_subs`,
	monic=`algcurves/monic`,
	almost_monic=`algcurves/amonic`,
	`puiseux/technical_answer`=`algcurves/puiseux_te`,
	`integral_basis/bound`=`algcurves/ib_bound`,
	Newtonpolygon=`algcurves/Newtonpolygon`
):

# This code uses the procedures in `algcurves/g_expand`.

macro(
	g_conversion1=`algcurves/g_conversion1`,
	g_conversion2=`algcurves/g_conversion2`,
	g_normal=`algcurves/g_normal`,
	g_expand=`algcurves/g_expand`,
	normal_tcoeff=`algcurves/normal_tcoeff`,
	g_factors=`algcurves/g_factors`,
	g_ext=`algcurves/g_ext`,
	truncate=`algcurves/truncate`
):

# This procedure lifts an expansion until it has multiplicity 1
# Input is a list v
# v[1] = The expansion so far
# v[2] = The accuracy, i.e. the lowest not yet determined coefficient. This
#  procedure starts with computing the coefficient of x^v[2]
# v[3] = Ramification, if v[2]=x^2, then x stands for sqrt(x).
# v[4] = The algebraic extensions in this expansion
# v[5] = The multiplicity of this factor
# v[6] = The degree of the algebraic extensions above the groundfield
# v[7] = The number Int_p, which is the sum of the valuations of the differences
#  of this expansion with the other expansions. This sum is needed for
#  computing an integral basis in an algebraic function field, cf JSC 94.
#  For an approximation r of a puiseux expansion p we have v(p-r)+Int_p=v(f(r))
# v[8] = list of the previous expansions, needed for determining the valuations
#   of the differences with the other Puiseux expansions
# f    = minimal polynomial of y over L(x), must be almost monic.
# Start computation with [0,0,x,ext,degree(f,y),1,1-degree(f,y),[]],f,x,y
lift_exp:=proc(v,f,x,y)
	local i,ii,r,res,v7,vv7,v3,ext,a,j,n,np,ram,j3;
	option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved. Author: M. van Hoeij`;
	if v[5]=1 then return {v} fi;
	v3:=degree(v[3],x); # the ramification index.
	res:={};
	r:=v[1]+y*x^v[2];
	vv7:=v[7]*v3+v[2]-1; # The ldegree of the previous step
	vv7:=vv7+v[5]; # Upperbound for the ldegree of this step
	# we try to find an equation for the unknown y
	ii:=truncate_subs(subs(x=v[3],f),x,y,r,vv7+1,v[4]);
	if ii=0 then
		error "degree estimate was wrong"
	fi;
	v7:=(ldegree(ii,x)-v[2])/v3; # The number Int_p=v(f(r))-v(p-r)
	r:=v_ext_m(g_factors(tcoeff(ii,x),y,v[4]),y);
	for i in r do
		res:=res union lift_exp([v[1]+x^v[2]*i[1],v[2]+1,v[3]
		 ,[op(i[3]),op(v[4])],i[2],v[6]*i[4],v7
		 ,[op(v[8]),[op(1..4,v)]]],f,x,y)
	od;
	if add(i[5]*i[6]/v[6]*degree(i[3],x)/v3,i=res)
	 <>degree(tcoeff(ii,x),y) then
		error "found wrong number of expansions"
	fi;
	if v[5]=degree(tcoeff(ii,x),y) then
		if ldegree(ii,x)<>vv7 then
			error "degree estimate was wrong"
		fi;
		return res
	fi;
	# Make sure that all valuations of the coefficients are correct
	ii:=collect(ii,y);
	ii:=add(normal_tcoeff(coeff(ii,y,i),x)*y^i
	 ,i=0..degree(ii,y));
	# Now compute all slopes s between 0 and 1
	np:=Newtonpolygon(ii,x,y);
	if nops(np)=2 and np[1][3]=0 then
		error "found wrong number of expansions"
	fi;
	for j in np do if nops(j)>2 and j[3]>0 and j[3]<1 then
		r:=g_factors(j[4],x,v[4]);
		r:=v_ext_m(r,x);
		for i in r do
			j3:=j[3]-v[2];
			ext:=[op(i[3]),op(v[4])];
			n:=mods(1/numer(j3),denom(j3));
			ram:=g_normal(i[1]^n)*x^denom(j3);
			a:=v[2]*denom(j3)-numer(j[3]);
			res:=res union lift_exp([
			 g_expand(subs(x=ram,v[1])+x^a*
			  g_normal(  i[1]^( (1-n*numer(j3))/denom(j3) )   )
			 ,ext),a+1,
			 g_expand(subs(x=ram,v[3]),ext),ext,i[2],i[4]*v[6],
			 (j[2]-j[1]*j[3]-a/degree(ram,x))/degree(v[3],x),
			 [op(v[8]),[op(1..4,v)]]],f,x,y)
		od
	fi od;
	res
end:

# This procedure lifts (= computes more terms of) expansions that have
# multiplicity 1. This could also be done with lift_exp, but this
# procedure is faster.
lift_exp_m1:=proc(v,f,x,y,nnk)
	local nk,res,powx,r,rr,v3,i;
	option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved. Author: M. van Hoeij`;
	v3:=degree(v[3],x); # The ramification index
	nk:=ceil(nnk*v3)/v3;
	powx:=v[2];
	if nk<powx/v3 then
		return [add(coeff(v[1],x,i)*x^i,
		 i=0..nk*v3-1),op(2..7,v)]
	elif nk=powx/v3 then
		return v
	fi;
	res:=v[1]+y*x^v[2];
	r:=truncate_subs(subs(x=v[3],f),x,y,res,v3*(nk+v[7]),v[4]);
	if ldegree(r,x)<>v[2]+v3*v[7] then
		error "wrong input"
	fi;
	r:=collect(r/x^ldegree(r,x),x);
	while nk>powx/v3 do
		rr:=collect(coeff(r,x,0),y);
		if degree(rr,y)<>1 then error "degree is not 1" fi;
		rr:=g_normal(-coeff(rr,y,0)/coeff(rr,y,1));
		res:=collect(subs(y=rr+x*y,res),x);
		powx:=powx+1;
		if nk<=powx/v3 then break fi;
		r:=subs(y=rr+x*y,r);
		r:=g_expand(r,v[4],x);
		if ldegree(r,x)<>1 then error "ldegree is not 1" fi;
		r:=add(coeff(r,x,i)*x^(i-1),i=1..v3*nk-powx)
	od;
	[subs(y=0,res),powx,op(3..7,v)]
end:

# Input: f in k[x,y]
# Output f(y=y_value) mod x^n
truncate_subs:=proc(f,x,y,y_value,n,ext)
	local ym,i,res,F;
	option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved. Author: M. van Hoeij`;
	if n=0 then
		return 0
	elif degree(f,x)>=n then
		F:=collect(f,x);
		F:=add(coeff(F,x,i)*x^i,i=0..n-1)
	else
		F:=f
	fi;
	res:=0;
	for i from 0 to degree(F,y) do
		if i=0 then ym:=1 else ym:=truncate(ym*y_value,n,x,ext) fi;
		res:=res+coeff(f,y,i)*ym
	od;
	res:=truncate(res,n,x,ext);
	# Make sure that ldegree(res,x) gives a right answer:
	normal_tcoeff(res,x)
end:

monic:=proc(f,x,y,ext,q)
	global g_conversion1,g_conversion2;
	local j,ff,lc,qq;
	option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved. Author: M. van Hoeij`;
	ff:=numer(g_normal(f));
	lc:=lcoeff(ff,y);
	if not has(lc,x) then
		q:=1;
		return g_normal(ff/lc)
	fi;
	lc:=subs(g_conversion1,sqrfree(subs(g_conversion2,lc),x));
	lc:=lc[2];
	lc:=g_expand(product(lc[j][1],j=1..nops(lc)),ext);
	ff:=monic(subs(y=y/lc,ff),x,y,ext,'qq');
	qq:=eval(qq);
	q:=g_expand(qq*lc,ext);
	ff
end:

# Make subs(x=0,lcoeff(f,y)) non-zero.
almost_monic:=proc(f,x,y,q)
	local ff,qq,i;
	option `Copyright (c) 1998 Waterloo Maple Inc. All rights reserved. Author: M. van Hoeij`;
	ff:=collect(evala(Primpart(f,y)),y);
	ff:=add(collect(coeff(ff,y,i),x,evala)*y^i,i=0..degree(ff,y));
	if coeff(lcoeff(ff,y),x,0)=0 then
		ff:=procname(subs(y=y/x,ff),x,y,'qq');
		q:=eval(qq)*x
	else
		q:=1
	fi;
	collect(ff,{x,y},distributed)
end:

# See ?puiseux
puiseux:=proc(aa,eqn::equation)
	global g_conversion1,g_conversion2;
	local x,v,a,f,y,ext,q,verz,i,res,n;
option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved. Author: M. van Hoeij`;
	if nargs<3 then error "too few arguments" fi;

if indets([args],RAD)<>{} then
	f:=[args];
	i := radfield(indets(f,RAD));
        return eval(subs(i[2],procname(op(eval(subs(i[1],f))))))
fi;
	if nargs=3 and not type(aa,RootOf) then
		error "wrong arguments"
	fi;
	x:=eval(op(1,eqn));
	v:=op(2,eqn);
	if v=infinity then
		a:=subs(x=1/x,aa)
	else
		a:=subs(x=x+v,aa)
	fi;
	ext:=g_ext([op(a),args[2..nargs]]);
if nargs=3 then
	f:=subs(g_conversion1,op(a));
	f:=subs(g_conversion2,subs(_Z=y,f));
	n:=args[3]
else
	f:=a;
	y:=args[3];
	n:=args[4]
fi;
	if member(n,{'minimal',`cycle structure`}) then n:=0 fi;
	if not type(n,numeric) then error "wrong arguments" fi;
	f:=almost_monic(f,x,y,'q');
	f:=subs(g_conversion1,f);
	if n<>0 then n:=n+ldegree(q,x) fi;
	verz:=`puiseux/technical_answer`(f,x,y,n,ext);
	if member(`cycle structure`,[args]) then
		return sort([seq(degree(i[3],x)$i[6],i=verz)])
	fi;
	res:={};
	for i in verz do
	  if nargs>4 then
		a:=subs(x=args[5],i[3]);
		res:=res union {[`if`(v=infinity, 1/x=1/a, x+v=a+v),
		 y=subs(x=args[5],i[1])/subs(x=a,q)]}
	  else
		res:=res union {subs(x=(x/lcoeff(i[3],x))^(1/degree(i[3],x)),
		 i[1])/q}
	  fi
	od;
	if v=infinity then
		res:=subs(x=1/x,res)
	else
		res:=subs(x=x-v,res)
	fi;
	subs(g_conversion2,res)
end:

# Note: the options remember won't work because of the local y in puiseux

# This procedure computes Puiseux expansions.
# It's input must be in internal format.
# Output: a list, see lift_exp
`puiseux/technical_answer`:=proc(f,x,y,n,ext)
	local i,v;
	options remember, 
	    `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved. Author: M. van Hoeij`;
	userinfo(5,'algcurves',`Computing Puiseux expansions...`);
	# It works for degree >=1, but not degree 0, so:
	if not has(f,y) then error "wrong input" fi;
if n=0 then
	i:=degree(collect(subs(x=0,f),y,g_normal),y);
	if i=-infinity then
		error "the curve should not have a factor involving only %1",x
	fi;
	i:=lift_exp([0,0,x,ext,i,1,1-i,[]],f,x,y)
else
	v:=procname(f,x,y,0,ext);
	if n=`integral basis` then
	  i:={seq(lift_exp_m1(i,f,x,y,`integral_basis/bound`(i,v,x)),i=v)}
	else
	  i:={seq(lift_exp_m1(i,f,x,y,n),i=v)}
	fi
fi;
	userinfo(5,'algcurves',`Done computing Puiseux expansions.`);
	i
end:

# Give an bound upper bound for the accuracy of the Puiseux
# expansion p that is required for the integral basis.
# Bound=max of (Int_i - Int_p + v(p-i)) taken over all Puiseux expansions i
#
# "An algorithm for computing an integral basis in an algebraic
#  function field", JSC 1994, 18, p. 353-363.
`integral_basis/bound`:=proc(p,v,x)
	local V,i,m,j,p1,p2,k;
	option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved. Author: M. van Hoeij`;
	# Note: in my paper it says: take the bound N_i > (Int_i - Int_p + v(p-i))
	# but the same proof shows that equality N_i = (Int_i - Int_p + v(p-i)) is
	# OK as long as we make sure that N_i > all v(p-i), which is done by the
	# following line:
	m:=p[2]/degree(p[3],x);
	p1:=[op(p[8]),p];
	for i in v minus {p} do
		# Now determine the valuation of the difference
		p2:=[op(i[8]),i];
		j:=1;
		while p1[j]=p2[j] do j:=j+1 od;
		V:=min(seq((k[2]-1)/degree(k[3],x),k=[p1[j],p2[j]]));
		m:=max(m,i[7]-p[7]+V)
	od;
	# Now m is the bound N_i for the expansion p.
	# A different bound is the highest multiplicity in the denominator
	# in the integral basis, which is given by the theorem in section 5
	# in my paper. Combine the bound N_i=m with this bound by taking the
	# minimum:
	min(m,floor(max(seq(i[7],i=v))))
end:

# Input: a monic polynomial f in k((x))[y]
# Output: sort of a Newton polygon
# Note: only the last element in the output is correct
Newtonpolygon:=proc(f,x,y)
	local i,n,ff,dummy,m,npg,powd,res,slope,val_powd,vals;
	options remember, 
	    `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved. Author: M. van Hoeij`;
	n:=degree(f,y);
	ff:=expand(f+(y+1)^n*x^(degree(f,x)+1));
	if n<2 then error "nothing to be factored" fi;
	vals:=[seq(ldegree(coeff(ff,y,i),x),i=0..n)];
	powd:=0;
	val_powd:=min(op(vals));
	npg:=[[powd,val_powd]];
	while powd<n do
		m:=min(seq((vals[powd+i+1]-val_powd)/i
		 ,i=1..n-powd));
		i:=n;
		while vals[i+1]-val_powd<>(i-powd)*m do i:=i-1 od;
		powd:=i;val_powd:=vals[i+1];
		npg:=[op(npg),[i,vals[i+1]]]
	od;
	res:=NULL;
	for i from 1 to nops(npg)-1 do
		# Now we compute the Newton polynomial of slope number i
		slope:=(npg[i+1][2]-npg[i][2])/(npg[i+1][1]-npg[i][1]);
		res:=res,[op(npg[i]),slope,add(
		 coeff(coeff(ff,y,denom(slope)*dummy+npg[i][1])
		  ,x,dummy*numer(slope)+npg[i][2])
		 *x^dummy,dummy=0..(npg[i+1][1]-npg[i][1])/denom(slope)
		)]
	od;
	[res,npg[nops(npg)]]
end:


#savelib ('lift_exp','lift_exp_m1','truncate_subs', \
 'almost_monic', \
 'monic','puiseux','`puiseux/technical_answer`', \
 '`integral_basis/bound`','Newtonpolygon'):