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

#--> differentials:  Compute a basis of the holomorphic differentials of an
#                    algebraic curve.
#
#  Input:        f   - A polynomial in 2 variables x and y, describing
#                      an irreducible algebraic curve C.
#               x,y  - variables
#
#  Calling sequence:
#		differentials(f,x,y);
#
#  Algorithm:   For n-multiple points you get n(n-1)/2 linear equations by
#		equating the ansatz A to zero at that point, and it's
#		derivatives of order < n-1.
#		For the remaining singularities, the special singularities,
#		where delta > n(n-1)/2, there are more linear conditions.
#		These are determined by a new method based on integral_basis.
#
#                                        Mark van Hoeij
#								June 1999.


macro(	differentials=`algcurves/differentials`,
	z=`algcurves/z`):

differentials:=proc(f,x::name,y::name)
	local a,d,ib,i,j,A,pr,r,eqn,ld,f_infty,AA,X,Y,k,var,W,Wi,S,v,g;
	global z,skip_dx;
	options remember,
`Copyright (c) 1999 Waterloo Maple Inc. All rights reserved. Author: M. van Hoeij`;
	if nargs<3 or not type(f,polynom(anything, [x, y])) then
		error "wrong number or type of arguments" 
	fi;
	g:=collect(primpart(f,{x,y}),{x,y},'distributed',Normalizer);
if g<>f then
	# Make sure input is normalized
	k:=procname(g,args[2..nargs])
elif not (nargs>3 and args[4]='skip_dx') then
	j:=d || x/diff(f,y);
	k:=[seq(i*j,i=procname(f,x,y,'skip_dx'))]
else
	d:=degree(f,{x,y});
	if (not has(f,x) or not has(f,y)) and d<>1 then
		error "%1 must have both variables %2, %3", f, x, y 
	elif d<3 then
		RETURN( [] )
	fi;
	userinfo(3,'algcurves',`Computing differentials...`);
	S:=algcurves['singularities'](f,x,y);
	g:=(d-1)*(d-2)/2-add(i[3]*`algcurves/degree_ext`(i,f),i=S);
	userinfo(1,'algcurves',`genus is`,g);
	userinfo(2,'algcurves',`Singularities with their mult and delta:`,S);
	v:={seq(`if`(i[3]>i[2]*(i[2]-1)/2,i[1],NULL),i=S)};
	userinfo(3,'algcurves',`Singularities where delta > mult*(mult-1)/2:`,v);
	if member([0, 1, 0],v) and member([1,0,0],v) then
		Wi:={[x,y],[y,x]}
	elif member([0, 1, 0],v) then
		Wi:={[y,x]}
	else
		Wi:={};
		for i in v do if i[3]=0 then Wi:={[x,y]} fi od
	fi;
	d:=d-3;
	k:=degree(f,x)-2;
	A:=add(add(a[i,j]*x^i*y^j,i=0..min(k,d-j)),j=0..min(d,degree(f,y)-2));
	var:=indets(A) minus {x,y};
	eqn:={};
	W:={};
	for i in S do
	if i[1]=[1,0,0] then
		AA:=expand(subs(x=1/z,y=y/z,A)*z^d);
		eqn:=eqn union
		 {seq(seq(coeff(coeff(AA,y,j),z,k-j),j=0..k),k=0..i[2]-2)}
	else
		if not member(i[1][3],{0,1}) or (i[1][3]=0 and i[1][2]<>1) then
			error "wrong syntax for singularities"
		fi;
		X,Y:=i[1][1],i[1][2];
		AA:=`if`(i[1][3]=0,expand(subs(x=x/z,y=1/z,A)*z^d),A);
		v:=map(evala@Expand,subs(x=X,y=Y,z=0,[AA,seq(seq(diff(
		 AA,`if`(i[1][3]=0,z,y)$j,x$(k-j)),j=0..k),k=1..i[2]-2)]));
		if i[3]<nops(v) then
			error "found wrong delta invariant"
		elif i[3]>nops(v) and i[1][3]=1 then
			W:={op(W),evala(Norm(x-X,{},indets(f,RootOf)))}
		fi;
		if type(Y,RootOf) and not member(Y,indets([f,X],RootOf)) then
			v:=map(coeffs,v,Y)
		fi;
		if type(X,RootOf) and not member(X,indets(f,RootOf)) then
			v:=map(coeffs,v,X)
		fi;
		eqn:=eqn union {op(v)}
	fi od;
	A:=subs(`solve/linear`(eqn,var),A);
	var:=indets(A) intersect var;
	eqn:={};
if W<>{} then
	# W contains x-coordinates of the finite special singularities:
	ib:=algcurves['integral_basis'](f,x,y,W);
	r:=evala(RootOf(f,y));
	for i in ib do
		pr:=subs(r=y,evala(Normal(subs(y=r,numer(i)*A),expanded)));
		pr:=evala(Expand(evala(Rem(numer(pr),denom(pr)*denom(i),x))));
		eqn:=eqn union map(coeffs,{coeffs(pr,y)},x)
	od
fi;
for k in Wi do
	# Determine equations for special singularities at infinity:
	X,Y:=op(k);
	f_infty:=subs(X=1,`algcurves/homogeneous`(f,X,Y,z,polynom));
	ib:=algcurves['integral_basis'](f_infty,z,Y,{[z,4]});
	AA:=expand(subs(X=1/z,Y=Y/z,A)*z^d);
	r:=evala(RootOf(f_infty,Y));
	for i in ib do
		pr:=subs(r=Y,evala(Normal(subs(Y=r,numer(i)*AA),expanded)));
		ld,pr:=ldegree(denom(pr)*denom(i),z),numer(pr);
		eqn:=eqn union
		 map(coeffs,{seq(expand(coeff(pr,z,j)),j=0..ld-1)},Y)
	od
od;
	A:=subs(`solve/linear`(eqn,var),A);
	var:=indets(A) intersect var;
	k:=solve(var);
	k:=[seq(subs(k,subs(i=1,A)),i=var)];
	if nops(k)<>g then
		WARNING("genus %1 <> nops(differentials), %2 must be reducible",g,f)
	fi;
	userinfo(3,'algcurves',`Done computing differentials`);
fi;
	# Try to simplify the output k, remove unnecessary factor -1.
	map(	proc(a) local b;
			b:=normal(a);b:=primpart(numer(b))/primpart(denom(b));
			`if`(op(1,b)=-1,-b,b)
		end
	,k)
end:

#savelib('differentials'):