# 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.
#
# This is a modification of the version found on:
#     http://www.math.fsu.edu/~hoeij/files/ConicProgram/
# The difference is that instead of trying one certificate, this
# version tries a number of certificates, computes a point for each,
# and then returns the smallest point that was found.
#
Conic := proc(a,b,c)
	local v,d,i,t,l,s,case,co,ansatz,c1,j,k,eq,f,n,si;

	# 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;
		v := map(ReduceList, [v]);
		v := sort(v,length)[1];

		# This type of sorting may not be necessary:
		#
		# v := map(i -> [op(i),length(factor(convert(i,`*`)))], v);
		# n := min(seq(i[4],i=v));
		# for i in v do if i[4]=n then v:=i; break fi od;

		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)];
	case := `if`(d mod 2 = [0,0,0], 0, 1);
	# The degrees for the ansatz:
	d := [seq(ceil(i/2)-case, 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 case = 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});
		f[k] := factors(primpart(v[k],t))[2];
		n[k] := nops(f[k])-`if`(case=0 or (k=3 and n[1]*n[2]>0), 0, 1);
		for l to nops(f[k]) do
			s := msqrt(v[i],v[j],f[k][l][1],t);
			if s = FAIL then return s fi;
			s := ansatz[i]*s[1]+ansatz[j]*s[2]*
				`if`(l>n[k],1,si[k,l]);
			eq := {op(eq),coeffs(collect(rem(s,f[k][l][1],t),t),t)}
		od
	od;

	d := {seq(seq(si[k,l],l=1..n[k]),k=1..3)};
	AllSols(eq, ansatz, indets(eq) minus indets([d,v]), d)
	fi
end:

AllSols := proc(eq, ansatz, var, si)
	local s,j;
	if si={} then
		# Solve linear equations, substitute the solution into ansatz.
		s := SolveTools:-Linear(eq, var);
		eval(ansatz, s)
	else
		s := si[1];
		seq(procname(eval(eq, s=j), ansatz, var, si minus {s}),j=[1,-1])
	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, x^2+factor(rem(a*b,p,t)));
	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), b*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: