# Input: a,b,c in Q(t1,t2,..,tk).
#
# Output: nonzero [x,y,z] for which a*x^2+b*y^2+c*z^2 = 0 if
#         such x,y,z exist in Q(t1,t2,..,tk) and FAIL otherwise.
#
# Author: Mark van Hoeij, Oct 2004.
#
Conic := proc(a,b,c)
	local v,d,i,t,l,s,d2,co,ansatz,c1,j,k,eq;

	# Some preprocessing:
	#
	v := [a,b,c];
	if   a=0 then [1,0,0]
	elif b=0 then [0,1,0]
	elif c=0 then [0,0,1]
	elif not type(v, list(ratpoly(rational))) then
		FAIL # Only Q(t1..tk) is implemented.
	elif args[-1]<>0 then
		d, v := ReduceConic(ReduceList(v));
		v := procname(op(v),0);
		if v = FAIL then return v fi;
		ReduceList([seq(v[i]/d[i], i=1..3)])
	elif indets(v)={} then
		# Solve a conic over Q (requires factoring integers)
		`algcurves/pyth_LLL`(op(v))
	else

	t := SelectVar(v);
	d := map(degree,v,t);
	l := map(lcoeff,v,t);
	s := d[1]+d[2]+d[3];
	d := [seq(s-d[i], i=1..3)];
	d2 := d mod 2;
	# The degrees for the ansatz:
	d := [seq(ceil(i/2)-max(op(d2)), i=d)];
	ansatz := [seq(add(co[i,j]*t^j, j=0..d[i]), i=1..3)];
	eq := {};

	# If the degrees are congruent mod 2, then we need some equations
	# coming from a conic with fewer variables:
	if d2 = [0,0,0] then
		s := procname(op(l));
		if s = FAIL then return s fi;
		# Equate leading coefficients of ansatz to entries of s.
		# Here c1 is a new variable to make the system homogeneous.
		eq := {seq(co[i,d[i]] - s[i]*c1, i=1..3)}
	fi;

	# Compute equations for each factor of a*b*c.
	for k to 3 do
		i,j := op({1,2,3} minus {k});
		for l in factors(primpart(v[k],t))[2] do
			s := msqrt(v[i],v[j],l[1],t);
			if s = FAIL then return s fi;
			s := ansatz[i]*s[1]+ansatz[j]*s[2];
			eq := {op(eq),coeffs(collect(rem(s,l[1],t),t),t)}
		od
	od;

	# Solve linear equations, substitute the solution into ansatz.
	s := SolveTools:-Linear(eq, indets(eq) minus indets(v));
	eval(ansatz, s);
	fi
end:


# Compute the square root of -a/b modulo p(t) and return
# a list with the numerator and denominator of this sqrt.
#
msqrt := proc(a,b,p,t)
	local x,R,v;
	if degree(b,t) > degree(a,t) then
		v := procname(b,a,p,t);
		return `if`(v=FAIL,v,[v[2],v[1]])
	fi;
	R := RootOf(p,t);
	v := subs(t = R, primpart(b*x^2+a,x));
	if degree(p,t)=1 then
		v := factors(v)
	else
		# If the number of vars is > 1 then the modular gcd code 
		# for function fields will really help speed things up here
		v := subs(R=t, evala(Factors(v, {R})))
	fi;
	v := select(has,v[2],x)[1][1];
	if degree(v,x)=2 then
		# Square root is non-rational:
		return FAIL
	fi;
	[coeff(v,x,0), coeff(v,x,1)]
end:

# Clear denominators and remove common factors of a list.
#
ReduceList := proc(v)
	local t,i,f;
	f := primpart(add(v[i]*t^i,i=1..nops(v)),t);
	[seq(coeff(f,t,i),i=1..nops(v))]
end:


# Select a good variable appearing in v.
#
SelectVar := proc(v)
	local vars,i,m,mn,mv;
	vars := indets(v);
	if vars = {} then error fi;
	for i in vars do
		m := max(op(map(degree,v,i)));
		if not assigned(mn) or m<mn then
			mn := m;
			mv := i
		fi;
	od;
	mv
end:

# Make a,b,c squarefree and coprime where v = [a,b,c].
#
ReduceConic := proc(v)
	local i,s,t,j,w,k,g;
	for i to 3 do
		s := sqrfree(v[i]);
		t[i] := mul(j[1]^iquo(j[2],2),j=s[2]);
		w[i] := mul(j[1]^irem(j[2],2),j=s[2])*s[1]
	od;
	for k to 3 do
		i,j := op({1,2,3} minus {k});
		g := gcd(w[i],w[j]);
		if g<>1 then
			w[i] := normal(w[i]/g);
			w[j] := normal(w[j]/g);
			w[k] := w[k]*g;
			t[k] := t[k]/g
		fi
	od;
	[t[1],t[2],t[3]], [w[1],w[2],w[3]]
end: