/* Perform a sequence of blow-ups on an algebraic curve to
make it non-singular */

define_variable(%npnum,0,any);

/* Calculate the vertices of the Newton polygon of the polynomial
p in variables x and y; return	a list of pairs [m[i],n[i]] of non-negative integers such that
	m[i]<m[i+1]
	n[i]>n[i+1]
	(m[i+1]-m[i])/(n[i+1]-n[i]) is increasing
and [m[i],n[i]] are the vertices of the Newton polygon spanned by L */

newton_polygon(p,x,y):=block([found,L],
	L:map(lambda([pow],[pow,lopow(coeff(p,x,pow),y)]),powers(p,x)),
	L:sort(L,lambda([p,q],p[1]<q[1])),
	while length(L)> 2 do (
	    for i:2 thru length(L)-1 do (
		found:false,
		if ( (l[i][2]-l[i-1][2])/(l[i][1]-l[i-1][1])>=
		   (l[i+1][2]-l[i][2])/(l[i+1][1]-l[i][1]) ) then (
			found:true,
			L:delete(L[i],L)
		)
	    ),
	    if not found then
		return()
	),
	if(length(L)=2 and L[1][2]=l[2][2])  then
		L:[first(L)],
	L
);

newton_powers(p,x,y):=block([L:newton_polygon(p,x,y)],
	makelist(-(l[i+1][2]-l[i][2])/(l[i+1][1]-l[i][1]),i,1,length(L)-1)
);

/* Return a list of substitutions which will blow up a polynomial
p in x and y at the point (0,0) */

np_substs(p,x,y):=block([L:newton_powers(p,x,y)],
	map(lambda([e],[x=x^num(e),y=x^denom(e)*y]),L)
);

/* Apply blowup substitutions to the polynomial p in x and y;
return a list of pairs [e,q] where e is a substitution and q is a
blown up polynomial p */

np_blowup(p,x,y):=
	map(lambda([e], block([q:subst(e,p)],[e,expand(q/x^lopow(q,x))])),
		np_substs(p,x,y));

/* Substitute x=0 and solve the resulting equation;
return a list [e,q], where e is a shift substitution and q
is the polynomial p with the substitution performed and expanded. */

np_shift(p,x,y):=
block([param,sol,s],
	%npnum:%npnum+1,
	param:concat("%np",%npnum),
	sol:ev(p/param^lopow(p,param),x=0,y=param,expand),
	sol:ev(sol/coeff(sol,param,hipow(sol,param)),expand),
	ev(tellrat(sol),algebraic:true),
	s:y=y+param,
	p:ev(rat(ev(p,s)),algebraic:true),
	[[sol=0,s],p]
);

np_resolution(p,x,y):=block([L:np_blowup(p,x,y)],
	display(L),
	L:makelist([e[1],np_shift(e[2],x,y)],e,L)
);