# Two conics over Q are isomorphic if there exists a birational map (an
# isomorphism) between them that is defined over Q. There is a simple way
# to test if two conics are isomorphic; just check if they have the
# same BadPrimes (also called ramified primes). The nontrivial part is to
# find an isomorphism.
#
# The main procedures in this file are:
#
# 1) Given two isomorphic conics, compute an isomorphism.
# 2) Given a conic, find a "simpler" conic that is isomorphic to it.
#
# See the file ExamplesIsom for the syntax of these procedures.
#
#            Author: Mark van Hoeij, Oct 2004


read pyth_LLL: # uncomment one line in that file if you have Maple 9

macro(
	SolveLinear = SolveTools:-Linear,
	IdentityMatrix = LinearAlgebra:-IdentityMatrix,
	parametrization = algcurves['parametrization']
):

#############################################################
#                 The main procedure                        #
#############################################################

ConicIsomorphism := proc(v,w)
	local V,W,mv,mw,P;
	if not assigned(_Env_ifactor_easy) then
		_Env_ifactor_easy := false
	fi;
	V := ReduceConic(op(v),mv);
	W := ReduceConic(op(w),mw);
	if assigned(_EnvConicConjugate) then
		# Set _EnvConicConjugate := [*,*,*] to get other isomorphisms
		mw := mw . SkewFieldConjugation(W)
	fi;
	if V=W then
		return mw . mv^(-1)
	fi;
	P := BadPrimes(V);
	if BadPrimes(W)<>P then
		error "No isomorphism exists"
	fi;
	if P=[] then
		# It's not really necessary to make a special case for this:
		P := P1ToConic(W) . P1ToConic(V)^(-1)
	else
		P := IsomReducedConics([V],[W],P)
	fi;
	P := mw . P . mv^(-1);
	if IsomorphismIsCorrect(P, v, w) then
		P
	else
		error "Found wrong answer"
	fi
end:

IsomorphismIsCorrect := proc(M, v, w)
	local fv,fw,x,y,z,P,i;
	fv := v[1]*x^2+v[2]*y^2+v[3]*z^2; # is 0 iff <x,y,z> is on conic v.
	P := M . <x,y,z>;  # Should be on conic w if <x,y,z> is on conic v.
	fw := expand(add(w[i]*P[i]^2,i=1..3));  # is 0 iff P is on conic w.
	testeq(rem(fw, fv, indets(fv)[1]))
end:

##################################################################
#                   Reduce a conic                               #
##################################################################

# Input: a,b,c which corresponds to the conic a*x^2+b*y^2+c*z^2 = 0.
# Output: the reduced conic, and a map that sends solutions of the
# new reduced conic to the old conic.
#
# Specifications:  First, remove denominators, then gcd's of coefficients,
# then remove square factors, then make all three positive or make
# two positive and one negative, then sort.
# If the optional argument M is given, then the transformation matrix
# will be assigned to it.
#
ReduceConic := proc(a::rational,b::rational,c::rational, M)
	options remember, system;
	local sa,sb,sc,A,B,C,g,j;
	if member(0,{a,b,c}) then
		error "Degenerate conic"
	fi;
	A := icontent(a+b*j+c*j^2);
	if convert(map(sign,[a,b,c]),`+`)<0 then A:=-A fi;
	A,B,C,sa,sb,sc := ConicReduceGCD(a/A, b/A, c/A);
	A,j := isq(A);
	sa := sa/j;
	B,j := isq(B);
	sb := sb/j;
	C,j := isq(C);
	sc := sc/j;
	g := <sa,0,0>,<0,sb,0>,<0,0,sc>;
	j := max(A,B,C);
	if A<>j then
		if B=j then
			A,B := B,A;
			g := g[2],g[1],g[3]
		else
			A,C := C,A;
			g := g[3],g[2],g[1]
		fi
	fi;
	if C>B then
		B,C := C,B;
		g := g[1],g[3],g[2]
	fi;
	HasBeenReduced(A,B,C) := true;
	if nargs=4 then M := Matrix([g]) fi;
	A,B,C
end:

HasBeenReduced := proc()
	options remember;
	# Can only reach this point if the input of this procedure
	# was not yet the output of ReduceConic.
	false
end:

ConicReduceGCD := proc(a,b,c)
	local A,B,C,sa,sb,sc,g1,g2,g3;
	A,B,C,sa,sb,sc := a,b,c,1,1,1;
	while {g1,g2,g3}<>{1} do
		g1 := igcd(A,B);
		A, B, C, sc := A/g1, B/g1, C*g1, sc*g1;
		g2 := igcd(B,C);
		A, B, C, sa := A*g2, B/g2, C/g2, sa*g2;
		g3 := igcd(C,A);
		A, B, C, sb := A/g3, B*g3, C/g3, sb*g3;
	od;
	A,B,C,sa,sb,sc
end:

##################################################################
#                         Notations                              #
##################################################################

# If the conic is [a,b,c] then the corresponding skew field F has
# basis b0,b1,b2,b3 where b0=1 and Q*b0 = Q is the center, b3 = c*b1*b2,
# bi*bj = -bj*bi for all 0<i<j, and b1^2 = -1/bc, b2^2 = -1/ac, b3^2 = -1/ab.
# Now a vector [x,y,z] corresponds to x*b1 + y*b2 + z*b3 in the
# skew field. Let v,w be two such vectors. Compute the product
# (without its b0-term) and return the corresponding vector

SkewFieldMultiplication := proc(a,b,c, v,w)
	[(v[2]*w[3]-v[3]*w[2])/a,  # b2*b3 = b1/a
	 (v[3]*w[1]-v[1]*w[3])/b,  # b3*b1 = b2/b
	 (v[1]*w[2]-v[2]*w[1])/c]  # b1*b2 = b3/c
end:

# The Norm (= -1 * the square) of [x,y,z] in the skew field F is
# (a*x^2+b*y^2+c*z^2)/abc  (times b0). Usually abc will be scaled away.

##################################################################
#          Find Isomorphism between reduced conics               #
##################################################################

# Compute the isomorphism when v,w are reduced conics, with bad primes P.
IsomReducedConics := proc(v,w,P)
	local V1,W1,V2,W2,N2,V3,W3,i;

	# Notation: Let Fv resp. Fw be the skew fields of v resp. w.
	# Take an element V1 in Fv and W1 in Fw with norm(V1)=norm(W1).
	#
	V1,W1 := FindMatchingNorm(v,w,P);

	# Take V2 in Fv and W2 in Fw with norm(V2)=norm(W2), such that
	# V2 resp. W2 anti-commutes with V1 resp W1. We just pick some W2
	# but then we need to compute to get a V2 with the same norm.
	#
	W2 := ChooseAntiCommute(W1);
	N2 := add(w[i]*W2[i]^2,i=1..3)/w[1]/w[2]/w[3]; # Norm of W2
	# Find a V2 that anti-commutes with V1, and has the same norm as W2:
	V2 := ComputeAntiCommute(op(v), V1, N2*v[1]*v[2]*v[3]);

	# Now multiply to find a third element of Fv resp. Fw
	#
	V3 := SkewFieldMultiplication(op(v), V1, V2);
	W3 := SkewFieldMultiplication(op(w), W1, W2);

	# Return a matrix that maps V1,V2,V3 to W1,W2,W3
	#
	Matrix([seq(<op(i)>,i=[W1,W2,W3])]) .
	Matrix([seq(<op(i)>,i=[V1,V2,V3])])^(-1)
end:
	
# Compute an element in skew fields of v resp. w that have the same norm.
FindMatchingNorm := proc(v,w,P)
	local N,W,i,j,vp,wp;
	N := {v[1]*v[2],v[1]*v[3],v[2]*v[3]} intersect
	     {w[1]*w[2],w[1]*w[3],w[2]*w[3]};
	if N<>{} then
		# We're in the easy case: One of the b.i's in Fv has
		# the same norm as one of the b.j's in Fw. All we
		# have to do is returning the corresponding vectors.
		N := N[1];
		vp := v[1]*v[2]*v[3];
		wp := w[1]*w[2]*w[3];
		for i to 3 do if vp = N*v[i] then break fi od;
		for j to 3 do if wp = N*w[j] then break fi od;
		vp := [0$3];
		subsop(i=1,vp), subsop(j=1,vp)
	else
		# The hard case. We just pick some W in Fw but then we'll
		# have to Find some element of Fv with the Same Norm:
		W,N := ChooseGoodElement(op(w), P);
		FindSameNorm(v,N*mul(v[i]/w[i],i=1..3),P), W
	fi:
end:

# v = [a,b,c]. Find [x,y,z] such that a*x^2+b*y^2+c*z^2 = N.
FindSameNorm := proc(v, N, Q)
	local P,e,inf,w,i,j;
	if not type(N,integer) then
		j := denom(N);
		return [seq(i/j, i=procname(v, N*j^2, Q))]
	fi;
	P := sort([op({op(Q),2,3,5,7,11,13,infinity})]);
	for e from 0 do
		inf := true;
		for i to 3 do
			w := subsop(i=e^2*v[i]-N, v);
			if BadPrimeInP(op(w),P,'inf') then
				next
			fi;
			w := pyth_LLL(op(w));
			if w<>FAIL then
				return [seq(`if`(j=i,e,w[j]/w[i]),j=1..3)]
			fi
		od;
		if inf then
			# Keep hitting bad prime infinity
			return [seq(i/2, i=procname(v, N*4, Q))]
		fi
	od
end:

# V is an element of the skew field. Find an element [x,y,z] that
# anti-commutes with V, and has norm N/abc, so ax^2+by^2+cz^2 = N.
ComputeAntiCommute := proc(a,b,c, V, N)
	local v,x,y,z,T,i;
	if length(V[1]) > min(length(V[2]),length(V[3])) or V[3] = 0 then
		v := procname(b,c,a,[V[2],V[3],V[1]], N);
		return [v[3],v[1],v[2]]
	fi;
	z := -(V[1]*a*x+V[2]*b*y)/V[3]/c;
	v := parametrization(expand(a*x^2+b*y^2+c*z^2-N),x,y,T);
	for i in [infinity,0,1,2,3]
	while not (type(x,rational) and type(y,rational)) do
		x := limit(v[1],T=i);
		y := limit(v[2],T=i)
	od;
	eval([x,y,z])
end:

##################################################################
#                 Choose simple vectors                          #
##################################################################

# Pick a small element in the skew field but do it
# in such a way that no odd bad primes appear in its norm.
# Output = small element of skew field, its norm
ChooseGoodElement := proc(a,b,c, P)
	local p;
	p := convert({op(P)} minus {2,infinity}, `*`);
	if igcd(c,p)=1 then
		[0,0,1], c
	elif igcd(b,p)=1 then
		[0,1,0], b
	elif igcd(a,p)=1 then
		[1,0,0], a
	else
		ifa(b+c); # store prime factors in b+c
		[0,1,1], b+c
	fi
end:

# W is an element of Fw chosen in ChooseGoodElement or in the easy case of
# FindMatchingNorm. Now choose an element that anti-commutes with W.
ChooseAntiCommute := proc(W)
	if W=[0,0,1] then
		[0,1,0]
	elif W[3]=0 then
		[0,0,1]
	elif W=[0,1,1] then
		[1,0,0]
	else
		error "Unexpected vector"
	fi
end:

##################################################################
#        Finding Bad Primes, and checking for Bad Primes         #
##################################################################

# The output is the list of all primes for which the conic has no solution
# over the completion of the rationals at that prime.
BadPrimes := proc(a,b,c)
	local v,w;
	options remember;
	if type(a,list) then
		return procname(op(a))
	elif member(0,{a,b,c}) then
		return []
	elif not HasBeenReduced(a,b,c) then
		return procname(ReduceConic(a,b,c))
	fi;
	w := seq(BadPrimesA(op(v)), v=[[a,b,c],[b,c,a],[c,a,b]]),
	 `if`(nops(map(sign,{a,b,c}))=1,infinity,NULL);
	if type(nops([w]),odd) then w := 2,w fi;
	sort([w])
end:

BadPrimesA := proc(a,b,c)
	local p,T;
	seq(`if`(Irreduc(b*T^2+c) mod p[1],p[1],NULL),p=ifa(a)[2])
end:

# Check if P contains a bad prime, where a,b,c need not be reduced.
# Stop as soon as a bad prime is encountered.
BadPrimeInP := proc(a,b,c,P,inf)
	local A,B,C,g,p,v,T,i,j;
	g := {a,b,c};
	if member(0,g) then
		error "input must be nonzero"
	elif member(infinity,P) and nops(map(sign,g))=1 then
		return true
	fi;
	if nargs=5 then inf := false fi;
	A,B,C := ConicReduceGCD(seq(numer(i)*denom(i), i=[a,b,c]))[1..3];
	if member(2,P) then
		A,B,C := seq(RemoveFactor(i,4),i=[A,B,C]);
		if irem(C,2)=0 then
			A,B,C := C,A,B
		fi;
		for i in [0,1,4] do for j in [0,1,4] do
			if irem(A*i+B*j+C,8)=0 then
				g := 0
			fi
		od od;
		if g<>0 then return true fi
	fi;
	for p in P do if p<>2 and p<>infinity then
		A,B,C := seq(RemoveFactor(i,p^2),i=[A,B,C]);
		v := subs(0=NULL,[A,B,C] mod p);
		if nops(v)=3 then next fi;
		if Irreduc(v[1]*T^2 + v[2]) mod p then
			return true
		fi
	fi od;
	false
end:

RemoveFactor := proc(a,f)
	local n;
	n:=a;
	while irem(n,f)=0 do n:=n/f od:
	n
end:

# Test if all primes in S are bad, where a,b,c are reduced.
# Stop as soon as a non-bad prime is encountered.
TestBadPrimes := proc(a,b,c,S)
	local x,y,z,f,p,ff;
	f := a*x^2+b*y^2+c*z^2;
	for p in S do
		ff := f mod p;
		ff := subs(indets(ff)[1]=1,ff);
		if not Irreduc(ff) mod p then
			return false
		fi
	od;
	true
end:

##################################################################
#           Finding a small equivalent conic                     #
##################################################################

# Find a simpler conic that is isomorphic (has the same BadPrimes). This 
# is not used by ConicIsomorphism, but can be combined nicely with it.
#
SimplifyConic := proc(a,b,c)
	local v,P,Q,i,S;
	options remember, system;
	if type(a,list) then
		return procname(op(a))
	elif not HasBeenReduced(a,b,c) then
		return procname(ReduceConic(a,b,c))
	fi;
	v := BadPrimes(a,b,c);
	if v=[] then return [1,1,-1] fi;
	S := {op(v)} minus {2,infinity};
	if member(infinity,v) then
		P := S
	else
		P := S union {-1}
	fi;
	for Q in [P, P union {2}] do for i in SortCombinations(Q) do
	if TestBadPrimes(op(i),S) then
		return i
	fi od od;
	# Found nothing simpler.
	sort([a,b,c],`>`)
end:

# Divide S into three pieces and take the product of each piece.
Combinations := proc(S)
	local e,i,r,m;
	if S={} then
		[1,1,1]
	elif nops(S)=1 then
		[S[1],1,1]
	elif nops(S)=2 then
		[S[1]*S[2],1,1], [S[1],S[2],1]
	else
		e := S[-1];
		if nops(S)>7 then
		  r := rand(1..3);
		  m := proc(L,e,a) subsop(a=e*L[a],L) end:
		  seq(m(i,e,r()),i=procname(S minus {e}))
		else
		  seq(op([[i[1]*e,i[2],i[3]],[i[1],e*i[2],i[3]],
		  [i[1],i[2],e*i[3]]]), i = procname(S minus {e}))
		fi
	fi
end:

SortCombinations := proc(Q)
	local v;
	v := [Combinations(Q)];
	v := map(sort,v,`>`);
	sort(v, proc(A,B)
		local a,b;
		a, b := member(-1,A), member(-1,B);
		if a=b then
		  a, b := member(1,A), member(1,B);
		  if a=b then
			return evalb(abs(A[2]+A[3])<abs(B[2]+B[3]))
		  fi
		fi;
		a
	end)
end:

##################################################################
#                    Auxiliary routines                          #
##################################################################

# Map the conic [1,1,-1] to the conic [a,b,c]. We could also have used
# IsomReducedConics instead, but this might be faster.
P1ToConic := proc(a,b,c)
	local v,x,y,t,R,i,j,a0,a1,a2,b0,b1,b2,c0,c1,c2;
	if [a,b,c]=[1,1,-1] then return IdentityMatrix(3) fi;
	v := parametrization(a*x^2+b*y^2+c,x,y,t);
	v := subs(t=(y+1)/x,v);
	v := primpart(v[1]+v[2]*t+t^2,t);
	R := (a0+a1*x+a2*y)*t^2 + (b0+b1*x+b2*y)*t + (c0+c1*x+c2*y);
	v := expand(rem(primpart(rem(v,R,t),t),x^2+y^2-1,y));
	v := SolveLinear({coeffs(v,{x,y,t})},indets(R) minus {t});
	R := expand(primpart(subs(v,R),t));
	R := [coeff(R,x,1), coeff(R,y,1), coeff(coeff(R,x,0),y,0)];
	Matrix([seq([seq(coeff(i,t,j),i=R)],j=0..2)])
end:

# Multiplying by this matrix replaces the three vectors W1,W2,W3 by
# conjugates in the skew field. One could take three variables c1,c2,c3
# and set _EnvConicConjugate := [c1,c2,c3]; and then afterwards figure out
# for which values of c1,c2,c3 the isomorphism matrix has small entries.
SkewFieldConjugation := proc(a,b,c)
	local x,y,z,d;
	x,y,z := op(_EnvConicConjugate);
	d := a*b*c+x^2*a+z^2*c+y^2*b;
	Matrix(3,3, [
	[(x^2*a-z^2*c+a*b*c-y^2*b)/d, -2*b*(-y*x+z*c)/d,  2*c*(z*x+y*b)/d],
	[2*a*(z*c+y*x)/d,   (a*b*c-x^2*a-z^2*c+y^2*b)/d, -2*c*(x*a-z*y)/d],
	[-2*a*(-z*x+y*b)/d, 2*b*(x*a+z*y)/d, (z^2*c-x^2*a+a*b*c-y^2*b)/d]])
end:

# This makes sure that algcurves[parametrization] uses pyth_LLL:
`algcurves/pyth2` := proc(f, x, y)
local P, c, i;
    if indets(f, {'name', 'RootOf'}) <> {x, y} or
    {degree(f, y), degree(f, {x, y})} <> {2} then error "wrong input"
    elif degree(f, x) = 1 then
        P := [subs(SolveLinear({subs(y = 0, f)}, {x}), x), 0]
    elif coeff(f, y, 1) <> 0 then
        c := -1/2*coeff(f, y, 1)/coeff(f, y, 2);
        P := procname(expand(subs(y = y + c, f)), x, y);
        if P <> FAIL then P := [P[1], P[2] + subs(x = P[1], c)] end if
    elif coeff(f, x, 1) <> 0 then
        c := -1/2*coeff(f, x, 1)/coeff(f, x, 2);
        P := procname(expand(subs(x = x + c, f)), x, y);
        if P <> FAIL then P := [P[1] + c, P[2]] end if
    else
        if not assigned(_Env_ifactor_easy) then _Env_ifactor_easy := true
        end if;
        P := primpart(`algcurves/homogeneous`(f, x, y, c));
        P := traperror(pyth_LLL(coeff(P, x, 2), coeff(P, y, 2),
            coeff(P, c, 2)));
        if P = lasterror or P = FAIL then P := FAIL
        else
            P := [P[1]/P[3], P[2]/P[3]];
            if subs(x = P[1], y = P[2], f) <> 0 then
                P := FAIL;
                userinfo(1, 'algcurves',
                    `Error in finding point on the conic`)
            end if
        end if
    end if;
    P
end proc: