##########################################################################
# This Maple routine computes sum of points on the jacobian of a
# genus 2 curve given by y^2=f[1]x^6+...+f[7]. Points on the jacobian
# are given by 4-tuples [x1,y1,x2,y2] where x1,y1 and x2,y2 are on
# the genus 2 curve. We may possibly have x=inf in which case we
# take for y the value y/x^3 rather than the infinite value of y.
# So on y^2=4x^6-1 the points at infinity are denoted by (inf,2) and 
# (inf,-2). On y^2=x^5-1 we have the single infinite point (inf,0).
#
# Usage: hypadd(f,p,q) where f=[f[1],...,f[7]] and p and q are 4-tuples
# representing pairs of points on the curve. The procedure returns
# another point on the jacobian in the form of a 4-tuple.
# The global variable inf is assumed unassigned.
#
# As an example consider the curve y^2=x^6-14x^4+49x^2-35. It contains
# for example the rational points (1,1),(1,-1),(2,1),(2,-1),(3,1),(3,-1)
# Some sample commands:
# hypadd([1,0,-14,0,49,0,-35],[1,1,1,1],[2,1,2,1]) yields [-3,1,-3,1]
# hypadd([1,0,-14,0,49,0,-35],[2,-1,2,-1],[3,1,inf,1]) 
# yields [213/28+sqrt(13393)/28, ...].
###########################################################################

hypadd:=proc(pol,p,q)
	local a,b,c,d,f,k,peq,\
	aa,bb,cc,x1,y1,x2,y2,x3,y3,x4,y4,\
	sol,spol,eq1,eq2,eq3,eq4,eta,eq,z,w,fpol;
	global x,inf;
	x:='x';
	peq:=a*x^3+b*x^2+c*x+d;
	f:=pol[1]*x^6+pol[2]*x^5+pol[3]*x^4+pol[4]*x^3+pol[5]*x^2+\
		pol[6]*x+pol[7];
	x1:=p[1];y1:=p[2];x2:=p[3];y2:=p[4];
	x3:=q[1];y3:=q[2];x4:=q[3];y4:=q[4];
#Check if points are really on the curve
	if not x1=inf and not expand(-y1^2+subs(x=x1,f))=0 then
		print(`x1 y1 not on curve`); RETURN() fi;
	if x1=inf and not expand(y1^2-pol[1])=0 then
		print(`x1 y1 not on curve`); RETURN() fi;
	if not x2=inf and not expand(-y2^2+subs(x=x2,f))=0 then
		print(`x2 y2 not on curve`); RETURN() fi;
	if x2=inf and not expand(y2^2-pol[1])=0 then
		print(`x2 y2 not on curve`); RETURN() fi;
	if not x3=inf and not expand(-y3^2+subs(x=x3,f))=0 then
		print(`x3 y3 not on curve`); RETURN() fi;
	if x3=inf and not expand(y3^2-pol[1])=0 then
		print(`x3 y3 not on curve`); RETURN() fi;
	if not x4=inf and not expand(-y4^2+subs(x=x4,f))=0 then
		print(`x4 y4 not on curve`); RETURN() fi;
	if x4=inf and not expand(y4^2-pol[1])=0 then
		print(`x4 y4 not on curve`); RETURN() fi;
	if y1=-y2 and x1=x2 then RETURN([x3,y3,x4,y4]) fi;
	if y1=-y3 and x1=x3 then RETURN([x2,y2,x4,y4]) fi;
	if y1=-y4 and x1=x4 then RETURN([x2,y2,x3,y3]) fi;
	if y2=-y3 and x2=x3 then RETURN([x1,y1,x4,y4]) fi;
	if y2=-y4 and x2=x4 then RETURN([x1,y1,x3,y3]) fi;
	if y3=-y4 and x3=x4 then RETURN([x1,y1,x2,y2]) fi;
#Now sort the points so that if [x1,y1],[x2,y2],[x3,y3],[x4,y4] has 
#repetitions then the distribution is either one of the following:
# (1234) (123,4) (13,24) (13,2,4) (12,34) (12,3,4) (1,2,3,4)
	if x1=x4 then x4:=x3;y4:=y3;x3:=x1;y3:=y1 fi;
	if x2=x3 then x2:=x1;y2:=y1;x1:=x3;y1:=y3 fi;
	if x2=x4 then x2:=x1;y2:=y1;x1:=x4;y1:=y4;\
		x4:=x3;y4:=y3;x3:=x1;y3:=y1 fi;
	if x3=x4 then x3:=x1;y3:=y1;x1:=x4;y1:=y4;\
		x4:=x2;y4:=y2;x2:=x1;y2:=y1 fi;
#Now we are going to make our equations
	if x1=x2 and x1=x3 and x1=x4 then
		eq:={op(makequation(f,peq,x1,y1,4))}
	fi;
	if x1=x2 and x1=x3 and not x1=x4 then 
		eq1:=makequation(f,peq,x4,y4,1);
		eq2:=makequation(f,peq,x1,y1,3);
		eq:={op(eq1),op(eq2)};
	fi;
	if x1=x3 and x2=x4 and not x1=x2 then
		eq1:=makequation(f,peq,x1,y1,2);
		eq2:=makequation(f,peq,x2,y2,2);
		eq:={eq1[1]+eq2[1],eq1[2]+eq2[2],\
			expand((x1-x2)*(eq1[1]-eq2[1])),\
			expand((x1-x2)*(eq1[2]-eq2[2]))}
	fi;
	if x1=x3 and not x1=x2 and not x2=x4 then
		eq1:=makequation(f,peq,x1,y1,2);
		eq2:=makequation(f,peq,x2,y2,1);
		eq3:=makequation(f,peq,x4,y4,1);
		eq:={op(eq1),op(eq2),op(eq3)}
	fi;
	if x1=x2 and x3=x4 and not x1=x3 then 
		eq1:=makequation(f,peq,x1,y1,2);
		eq2:=makequation(f,peq,x3,y3,2);
		eq:={op(eq1),op(eq2)}
	fi;
	if x1=x2 and not x3=x4 and not x1=x3 then 
		eq1:=makequation(f,peq,x1,y1,2);
		eq2:=makequation(f,peq,x3,y3,1);
		eq3:=makequation(f,peq,x4,y4,1);
		eq:={op(eq1),eq2[1]+eq3[1],\
			expand((x3-x4)*(eq2[1]-eq3[1]))}
	fi;
	if not x1=x3 and not x1=x2 then 
		eq1:=makequation(f,peq,x1,y1,1);
		eq2:=makequation(f,peq,x2,y2,1);
		eq3:=makequation(f,peq,x3,y3,1);
		eq4:=makequation(f,peq,x4,y4,1);
		eq:={eq1[1]+eq2[1],expand((x1-x2)*(eq1[1]-eq2[1])),\
			eq3[1]+eq4[1],expand((x3-x4)*(eq3[1]-eq4[1]))}
	fi;
	sol:=op(solve(eq,{a,b,c,d}));
	spol:=subs(sol,peq);
	fpol:=expand(f-spol^2);
	if not x1=inf and not x2=inf then
		eta:=expand((x-x1)*(x-x2));
		fpol:=simplify(fpol/eta)
	fi;
	if x1=inf and not x2=inf then
		fpol:=simplify(fpol/(x-x2))
	fi;
	if not x1=inf and x2=inf then
		fpol:=simplify(fpol/(x-x1))
	fi;
	if not x3=inf and not x4=inf then
		eta:=expand((x-x3)*(x-x4));
		fpol:=simplify(fpol/eta)
	fi;
	if x3=inf and not x4=inf then
		fpol:=simplify(fpol/(x-x4))
	fi;
	if not x3=inf and x4=inf then
		fpol:=simplify(fpol/(x-x3))
	fi;
	if coeff(fpol,x,2)=0 and coeff(fpol,x,1)=0 then
		RETURN([inf,-coeff(spol,x,3),inf,-coeff(spol,x,3)]);
	elif coeff(fpol,x,2)=0 then
		z:=-coeff(fpol,x,0)/coeff(fpol,x,1);
		w:=subs(x=z,spol);
		RETURN([inf,-coeff(spol,x,3),z,-w]);
	else
		z:=[solve(fpol,x)];
		w:=[expand(subs(x=z[1],spol)),expand(subs(x=z[2],spol))];
		RETURN([z[1],-w[1],z[2],-w[2]])
	fi end:

makequation:=proc(g,p,xi,eta,n)
	local h,q,t,l,h1,h2,h3,q0,q1,q2,q3,k;
	global x;
	if xi=inf then
		h:=expand(t^6*subs(x=1/t,g));
		q:=expand(t^3*subs(x=1/t,p));
	else
		h:=expand(subs(x=t+xi,g));
		q:=expand(subs(x=t+xi,p));
	fi;
	h1:=coeff(h,t,1);
	h2:=coeff(h,t,2);
	h3:=coeff(h,t,3);
	q0:=coeff(q,t,0);
	q1:=coeff(q,t,1);
	q2:=coeff(q,t,2);
	q3:=coeff(q,t,3);
	l:=[q0-eta];
	if n>1 then l:=[op(l),2*eta*q1-h1] fi;
	if n>2 then l:=[op(l),(2*eta)^3*q2-(2*eta)^2*h2+h1^2] fi;
	if n>3 then l:=[op(l),\
		16*eta^5*q3-8*eta^4*h3+4*eta^2*h1*h2-h1^3] fi;
	RETURN(l) end: