# Input: x = [x1, x2, x3]
# q is an irreducible homogeneous of degree 2 in x1, x2, x3
# p = [p1, p2, p3] is a point on q.
#
# The output is a list X such that the substitution {x[1]=X[1], x[2]=X[2], x[3]=X[3]}
# turns q to the form x2^2 - x1*x3
#
Simplify_q := proc(q, x::list, p::list, X::list)
	local x1,x2,x3,Q,i,pts,Q1,Q2,Q3,Line,Lines,EQ,v,SL,S;
	if nops(x) <> 3 or nops(p) <> 3 or p = [0,0,0] then
		error "wrong input"
	elif nargs = 3 then
		return procname(args, x)
	fi;
	x1,x2,x3 := op(x);
	if p[3]=0 then
		return procname(q, x, [p[3],p[1],p[2]], subs({x1=x2,x2=x3,x3=x1},X))
	elif p[3] <> 1 then
		return procname(q, x, normal([p[1]/p[3], p[2]/p[3], 1]), X)
	elif p[1] <> 0 or p[2] <> 0 then
		return procname(q, x, [0, 0, 1], subs(x1=x1 + p[1]*x3, x2=x2 + p[2]*x3, X))
	fi;
	Q := collect(subs({x1=X[1],x2=X[2],x3=X[3]}, q), {x1,x2,x3}, normal);
	if normal({seq(coeff(Q,x[i],2),i=1..3)}) <> {0} then
		pts := {[0,0,1]};
		# Now Q has the point [0,0,1] on it.  Lets intersect Q with three lines through
		# that point to find more points, until we have three points.
		Q1 := normal(subs(x1 = 0, Q)/x2);
		if coeff(Q1,x3)<>0 then
			pts := pts union {[0, 1, normal(-coeff(Q1,x2)/coeff(Q1,x3))]}
		fi;
		Q2 := normal(subs(x2 = 0, Q)/x1);
		if coeff(Q2,x3)<>0 then
			pts := pts union {[1, 0, normal(-coeff(Q2,x1)/coeff(Q2,x3))]}
		fi;
		Q3 := normal(subs(x1 = x2, Q)/x2);
		if nops(pts) < 3 and coeff(Q3,x3)<>0 then
			pts := pts union {[1, 1, normal(-coeff(Q3,x2)/coeff(Q3,x3))]}
		fi;
		if nops(pts) < 3 then
			error "Input must be irreducible, not a product of two lines"
		fi;
		Line := a1 * x1 + a2 * x2 + a3 * x3;
		Lines := {};
		for i in pts do
			EQ := {seq(eval(Line, {x1 = v[1], x2 = v[2], x3 = v[3]}), v = pts minus {i})};
			Lines := Lines union {primpart(subs(solve(EQ, {a1,a2,a3}), Line), {x1,x2,x3})}
		od;
		SL := solve({seq(y[i] = Lines[i], i=1..3)}, {x1,x2,x3});
		return procname(q, x, p, subs(seq(y[i]=x[i],i=1..3), subs(SL, X)))
	fi;
	# Now [1,0,0], [0,1,0], [0,0,1] are on the curve.
	# So we're looking at .. x1x2 + ..x1x3 + ..x2x3 and all of those coefficients
	# are non-zero (otherwise the equation would be reducible).
	# Aiming for x2^2 - x1 * x3
	S := subs(x1 = x1 - x2 * coeff(coeff(Q,x2),x3) / coeff(coeff(Q,x1),x3), X);
	Q := collect(subs({x1=S[1],x2=S[2],x3=S[3]}, q), {x1,x2,x3}, normal);
	S := subs(x3 = -x3*lcoeff(Q,x2)/coeff(coeff(Q,x3),x1), S);
	Q := collect(subs({x1=S[1],x2=S[2],x3=S[3]}, q), {x1,x2,x3}, normal);
	S := subs(x3 = x3 + x2 * coeff(coeff(Q,x2),x1)/lcoeff(Q,x2), S);
	Q := primpart(collect(subs({x1=S[1],x2=S[2],x3=S[3]}, q), {x1,x2,x3}, normal), x2);
	if not member(x2^2-x1*x3, {Q, -Q}) then error "unexpected" fi;
	S;
end:

q := 3*x1^2+x1*x2+x1*x3+x2^2-10*x2*x3-20*x3^2;
Simplify_q(q, [x1,x2,x3], [3,2,1]);

primpart(subs( {x1 = %[1], x2=%[2], x3=%[3]}, q), x2);