#Implementation of some algorithms by that new guy, C F Gauss.

#His 'reduction' amounts to finding a matrix  C  with  C Q C^t  being the
#symmetric matrix representing his 'reduced' quadratic form.
#All our reduction procedures will return the corresponding  C  rather
#than simply the reduced matrix  C Q C^t, which is then easily found anyway.
#So general procedure is: given Q, find C1, make Q1, induct C0, return C0 C1.

#He then applies this to solving Legendre's equation: he calculates a change
#of basis in this way so that the quadratic form  ax^2+by^2+cz^2 become just
#(-abc)(X^2+YZ), two solutions of which are clearly [X,Y,Z]=[0,1,0] and [0,0,1]
#which are then translated back to xyz-coordinates.

#How good is this solution? Ingenious, absolutely, but terrible in practice.
#Two impediments:
#1. 3D-lattice reduction requires reducing 2D lattices with _indefinite_ forms.
#  Gauss's algorithm is NOT polynomial-time; I see no way around that, either.
#2. An initial call to "Lemma279" greatly increases the magnitude of the
#  numbers involved. Probably this part can be fixed without much ingenuity.

#Warning: his matrices "S" seem to be my  C^t, not C. And I may not have his
#definition of 'reduced' quite matched. If I were really to
#use this stuff, I would return [C,CQC^t] instead of just C. Also I'd
#switch to matrix input rather than list input. I think.

##############################################################################

with(linalg):
I2:=matrix(2,2,[1,0,0,1]):
I3:=matrix(3,3,[1,0,0,0,1,0,0,0,1]):
J2:=matrix(2,2,[0,1,1,0]):

#Maple's isolve() does not return the same answer all the time,
# thus affecting predictability of answers! So instead I use
#
euc2:=proc(a,b) local k,x:
if abs(a)<abs(b) then x:=euc2(b,a): RETURN([x[2],x[1]]): fi:
if b=0 then RETURN([sign(a),0,0]):fi:
k:=round(a/b): x:=euc2(a-k*b,b): RETURN([x[1],x[2]-k*x[1]]): end:

chrem2:=proc(a,b)
# Maple idiocy: amazingly, chrem([a1,a2],[1,c]) is OK but not the other way!
if abs(b[2])=1 then (a[1] mod b[1]) else chrem(a,b) fi;
end:

################################

#Binary forms stuff first. 
#Several ways (mind,minn,minq,minz) to reduce quadratic forms.
#Gauss writes (a,b,c) for the quadratic form with matrix Q=[[a,b],[b,c]].

#For definite forms, this is ordinary lattice reduction: (or, see Art. 171)
#Here I imported an old routine so its recursion is a little different.
#
mind:=proc(v1,v2,Q1)    
#Assume  v1, v2 lie in  Z^2  and  Q1=[a,b,c] the Gauss notation for a definite
#quadratic form on Z^2: Q1(x,y)=ax^2+2bxy+cy^2 (art 153).
#Return coordinates [x,y,z,w] of a good basis  {x v1 + y v2, z v1 + w v2}.
local L1,L2,IP,sol,k,Q,M1:
#Q must be _positive_ definite:
if Q1[1]<0 then Q:=[-Q1[1],-Q1[2],-Q1[3]] else Q:=Q1 fi:
L1:=Q[1]*v1[1]^2+2*Q[2]*v1[1]*v1[2]+Q[3]*v1[2]^2:
L2:=Q[1]*v2[1]^2+2*Q[2]*v2[1]*v2[2]+Q[3]*v2[2]^2:
if L2<L1 then sol:=mind(v2,v1,Q):RETURN(multiply(sol,J2)) fi:
IP:=Q[1]*v1[1]*v2[1]+Q[2]*(v1[1]*v2[2]+v2[1]*v1[2])+Q[3]*v1[2]*v2[2]:
k:=round(IP/L1):
if k=0 then RETURN(I2):fi:
if k=1 and IP/L1=1/2 then RETURN(I2):fi:
sol:=mind(expand(v2-k*v1),v1,Q):M1:=matrix(2,2,[-k,1,1,0]):
RETURN(multiply(sol,M1)):
end:

#For non-definite forms, this is a non-unique reduction! See art 183 et seq
#Warning: this implementation (Gauss's own) is NOT effective.
#
minn:=proc(a,b,c)
# Reduce indefinite form so | |A|-sqrt(d) |<B < sqrt(d), d=B^2-AC. 
local a1,b1,a2,d,k,M1,M2:
lprint(`In minn with`,a):
a1:=c:d:=b^2-a*a1:
#square discriminants are different:
if isqrt(d)^2=d then RETURN(minq(a,b,c)): fi:
if b^2<d and b>0 and b<abs(a) and abs(d-a^2-b^2)<2*abs(a)*b
  then RETURN(I2): fi:
b1:=(-b) mod a1:
k:=round(((isqrt(d)-b1)/abs(a1)-1/2)):if abs(a1)*(k+1)+b1<0 or (abs(a1)*(k+1)+b1)^2 < d then k:=k+1: print(`rounding fix`):fi:
b1:=b1+abs(a1)*k:
a2:=(b1^2-d)/a1:
if abs(a2)<abs(a1) then M1:=minn(a1,b1,a2): else M1:=I2: fi:
k:=-(b+b1)/a1;
M2:=matrix(2,2,[0,-1,1,k]):
RETURN(multiply(M1,M2)):
end:
#Examples with many steps arise e.g. cc=rand(): ifactor(cc^2-29):aa:=...
#close to sqrt; cc:=""/": mywin(aa,bb,cc):
#24669869   -23575126   22528963

#Forms with a square "determinant" are in art 206:
#
minq:=proc(a,b,c)
local h,alp,bet,gam,del,su,a2,A,k,ans,M1,M2:
if a=0 and c=0 then RETURN(I2): fi:
if a=0 then ans:=minq(c,b,a):RETURN(multiply(ans,J2)) fi:
h:=sqrt(b^2-a*c):
del:=denom((h-b)/a):bet:=del*(h-b)/a:
euc2(del,-bet):alp:="[1]:gam:=""[2]:
M1:=matrix(2,2,[alp,gam,bet,del]):
#What if not yet good enough? transform more; another matrix to the left
a2:=a*alp^2+2*b*alp*gam+c*gam^2:A:=(a2 mod 2*h):k:=(A-a2)/2/h:
M2:=matrix(2,2,[1,k,0,1]):
RETURN(multiply(M2,M1)):
end:

#For forms with discriminant zero see article 215 
#
minz:=proc(a,b,c)  
local m,g,h,u,v:
m:=gcd(a,c):if a<0 then m:=-m fi:if c*m<0 then m:=-m fi:
g:=sqrt(a/m):h:=sqrt(c/m):if b*m < 0 then h:=-h:fi:
euc2(g,h):u:="[1]:v:=""[2]:
RETURN(matrix(2,2,[h,-g,u,v])):
end:

#So here we just say "reduce it!" for any binary quadratic form
#
binmin:=proc(a,b,c) local d,ans,X,Y,x,y:
d:=b^2-a*c:
if d<0 then ans:=mind([1,0],[0,1],[a,b,c]) 
	elif d>0 then ans:=minn(a,b,c):
	else ans:=minz(a,b,c) fi:
RETURN(ans):
end:

#A checker:
#
checkbin:=proc(C,Q) print(multiply(multiply(C,matrix(2,2,[Q[1],Q[2],Q[2],Q[3]])),transpose(C)));end: #should show nice,small symmetric matrix
#
#Or just see the reduced form:
red:=proc(a,b,c) checkbin(binmin(a,b,c),[a,b,c]): end:

################################

#We follow Gauss' discussion of ternary quadratic forms, Art 266 et seq.
#He represents the quadratic form given by a symmetric matrix  Q  by
#the sextuple  [Q11,Q22,Q33,Q23,Q13,Q12]. Note his 'determinant' is -det(Q).

#art 267:
adjo:=proc(a,a1,a2,b,b1,b2)
[b^2-a1*a2,b1^2-a*a2,b2^2-a*a1,a*b-b1*b2,a1*b1-b*b2,a2*b2-b*b1]
end:
#The adjoint matrix of  Q  is again symmetric

#Now we consider reductions of ternary forms
#
#art 272: He proposes to alternate the reduction of binary quadratic forms
#on the subspaces  span(e1,e2) and span(e2,e3). 

min1:=proc(a,b,c,d,e,f) local ans,a1,b1,c1,d1,e1,f1,M1,M2,Q,Q2:
lprint(`1:`,a,b,c,d,e,f,a*b*c+2*d*e*f-a*d^2-b*e^2-c*f^2);
ans:=binmin(a,f,b): 
if ans=I2 then RETURN(I3) fi:
M1:=matrix(3,3,[ans[1,1],ans[1,2],0,ans[2,1],ans[2,2],0,0,0,1]):
Q:=matrix(3,3,[a,f,e,f,b,d,e,d,c]):
Q2:=multiply(multiply(M1,Q),transpose(M1)):
a1:=Q2[1,1]:b1:=Q2[2,2]:c1:=Q2[3,3]:d1:=Q2[2,3]:e1:=Q2[3,1]:f1:=Q2[1,2]:
if [a,b,c,d,e,f]=[a1,b1,c1,d1,e1,f1] then RETURN(I3) fi:
M2:=min2(a1,b1,c1,d1,e1,f1):
RETURN(multiply(M2,M1)):
end:

min2:=proc(a,b,c,d,e,f) local ans,a1,b1,c1,d1,e1,f1,M1,M2,Q,Q2,d0:
lprint(`2:`,a,b,c,d,e,f,a*b*c+2*d*e*f-a*d^2-b*e^2-c*f^2);
ans:=binmin(f^2-a*b,a*d-e*f,e^2-a*c):
if ans=I2 then RETURN(I3) fi:
d0:=ans[1,1]*ans[2,2]-ans[1,2]*ans[2,1]:#=+-1
M1:=matrix(3,3,[d0,0,0,0,ans[1,1],-ans[1,2],0,-ans[2,1],ans[2,2]]):#inverse transf.
Q:=matrix(3,3,[a,f,e,f,b,d,e,d,c]):
Q2:=multiply(multiply(M1,Q),transpose(M1)):
a1:=Q2[1,1]:b1:=Q2[2,2]:c1:=Q2[3,3]:d1:=Q2[2,3]:e1:=Q2[3,1]:f1:=Q2[1,2]:
M2:=matrix(3,3,min1(a1,b1,c1,d1,e1,f1)):
RETURN(multiply(M2,M1)):
end:

#Further reductions from Art 274.
#
minx:=proc(a,b,c,d,e,f)
local r,s,t,A,k,n,m,N:
lprint(`3:`,a,b,c,d,e,f):
A:=adjo(a,b,c,d,e,f):
if a*A[3]=0 and a+A[3]<>0 then 
 print(`See Art. 274: have a A2=0 but not both zero`):
 RETURN([0,0,0]):
fi:
if a=0 and A[3]=0 then
 k:=gcd(b,e):n:=d mod k: if n>k/2 then n:=n-k: fi:
 euc2(b,e):t:="[1]:r:=""[2]:t:=t*(n-d)/k:r:=r*(n-d)/k:
 m:=(c+2*d*t+a*t^2) mod 2*e:
 if m>abs(e) then m:=m-2*abs(e) fi:
 s:=(m-(c+2*d*t+a*t^2))/2/e:
else
 r:=round(-f/a):
 t:=round(A[4]/A[3]):N:=A[4]-A[3]*t:
 s:=round((A[5]-N*r)/A[3]):
fi:
RETURN([r,s,t]):
end:

#We polish off unimodular forms with Art. 277
#
minu:=proc(a,b,c,d,e,f)
local dd:
if d=0 and e=0 and f=0 then #form is diagonal; switch to +-(x^2+2yz) if indef.
 if a=b and b=c then RETURN(I3) fi: #that's for definite unimodular forms.
 dd:=-a*b*c:
 if a=b then RETURN(matrix(3,3,[1,1,1,dd,0,dd,0,-1,-1]));
	elif b=c then RETURN(matrix(3,3,[1,1,1,dd,0,dd,-1,-1,0])) 
	else RETURN(matrix(3,3,[1,1,1,dd,dd,0,0,-1,-1])):fi:
fi:
#Next are the cases of column-permutations of Gauss's "(+-1, 0, 0; +-1, 0, 0)"
if [b,c,e,f]=[0,0,0,0] then RETURN(matrix(3,3,[1,0,0,0,1,0,0,0,d])) fi:
if [a,c,d,f]=[0,0,0,0] then RETURN(matrix(3,3,[0,1,0,0,0,1,e,0,0])) fi:
if [a,b,d,e]=[0,0,0,0] then RETURN(matrix(3,3,[0,0,1,1,0,0,0,f,0])) fi:
#Now the "(0,+-1,+-1; 0,+-1,0)" 's
if [a,d,f]=[0,0,0] then RETURN(matrix(3,3,[-c*e,1,0,-e*(c+b)/2,b*c,1,e,0,0])):fi:
if [b,d,e]=[0,0,0] then RETURN(matrix(3,3,[0,-a*f,1,1,-f*(a+c)/2,c*a,0,f,0])):fi:
if [c,e,f]=[0,0,0] then RETURN(matrix(3,3,[1,0,-b*d,a*b,1,-d*(b+a)/2,0,0,d])):fi:
lprint(a,b,c,d,e,f,`unimodular reduction failure!`):
end:

#So we just call this:
#
ternmin:=proc(a,b,c,d,e,f) local M1,M2,M3,Q,Q2,Q3,a1,b1,c1,d1,e1,f1:
M1:=min1(a,b,c,d,e,f): if M1=I3 then M1:=min2(a,b,c,d,e,f): fi:
Q:=matrix(3,3,[a,f,e,f,b,d,e,d,c]):
Q2:=multiply(multiply(M1,Q),transpose(M1)):
a1:=Q2[1,1]:b1:=Q2[2,2]:c1:=Q2[3,3]:d1:=Q2[2,3]:e1:=Q2[3,1]:f1:=Q2[1,2]:
minx(a1,b1,c1,d1,e1,f1):
M2:=matrix(3,3,[1,0,0,"[1],1,0,"[2],"[3],1]):
if abs(a*b*c+2*d*e*f-a*d^2-b*e^2-c*f^2)=1 then
   Q3:=multiply(multiply(M2,Q2),transpose(M2)):
   a1:=Q3[1,1]:b1:=Q3[2,2]:c1:=Q3[3,3]:d1:=Q3[2,3]:e1:=Q3[3,1]:f1:=Q3[1,2]:
   M3:=minu(a1,b1,c1,d1,e1,f1):
   RETURN(multiply(M3,multiply(M2,M1))):
else RETURN(multiply(M2,M1)): fi:
end:

#A double-check:
#
#with(linalg):
checktern:=proc(CC,Q) local AA,BB:
BB:=matrix(3,3,[Q[1],Q[6],Q[5],Q[6],Q[2],Q[4],Q[5],Q[4],Q[3]]):
print(det(CC));  #get 1 (or -1?)
AA:=multiply(multiply(CC,BB),transpose(CC)); 
print(AA):       #get simple symmetric matrix
[AA[1,1],AA[2,2],AA[3,3],AA[2,3],AA[3,1],AA[1,2]]:
end:
#usage:  checktern(",Q);   right after a ternmin().

################################

#Now we solve Legendre's equation

#First some prep: (Art. 279)

euc3:=proc(a,b,c) local k,x:
if abs(a)<abs(b) then x:=euc3(b,a,c): RETURN([x[2],x[1],x[3]]): fi:
if abs(b)<abs(c) then x:=euc3(a,c,b): RETURN([x[1],x[3],x[2]]): fi:
if b=0 then RETURN([sign(a),0,0]):fi:
k:=round(a/b): x:=euc3(a-k*b,b,c): RETURN([x[1],x[2]-k*x[1],x[3]]): end:

l279:=proc(a1,a2,a3) #complete [a1,a2,a3,???,???,...] to a unimodular matrix
# Warning: this implementation (Gauss's) produces large numbers; it is not
#  clear how easily this can be fixed -- how to choose better  e1,e2,e3?
local A1,A2,A3,B1,B2,B3,C1,C2,C3,alp,bet,As,es,e1,e2,e3,b1,b2,b3,h:
As:=euc3(a1,a2,a3); alp:=As[1]*a1+As[2]*a2+As[3]*a3:
A1:=As[1]:A2:=As[2]:A3:=As[3]:
e1:=1:e2:=0:e3:=0: if A2=0 and A3=0 then e2:=1 fi:
b1:=e2*A3-e3*A2:b2:=e3*A1-e1*A3:b3:=e1*A2-e2*A1:
bet:=gcd(gcd(b1,b2),b3);
C1:=(a2*b3-a3*b2)/alp/bet;
C2:=(a3*b1-a1*b3)/alp/bet;
C3:=(a1*b2-a2*b1)/alp/bet;
es:=euc3(b1,b2,b3):
e1:=es[1]:e2:=es[2]:e3:=es[3]:
h:=e1*a1+e2*a2+e3*a3;
B1:=alp*e1-h*A1;B2:=alp*e2-h*A2;B3:=alp*e3-h*A3;
[B1,B2,B3,C1,C2,C3];
end:

#Legendre's theorem, Art 294 (First proof; second, art 295, not understood yet)

gsolve:=proc(a,b,c,nn,p,q)
# The point of this proof is to use  nn,p,q  to find a change-of-coordinate
# matrix which transforms the quad.form ax^2+by^2+cz^2 to (abc)(x^2+yz); then
# easy solutions arise from taking xyz=[0,1,0] or [0,0,1], i.e. the last two
# rows of the output.
local aa,bb,cc,A,B,C,mu,As,A1s,g1,d,M3,Q,M2:
d:=-a*b*c;
# First reduce to an equivalent quadratic problem with a quadratic form g which
# includes summands ...+m'y^2+2nyz+m''z^2+... with d | {m', n, m''}
Q:=matrix(3,3,[a,0,0,0,b,0,0,0,c]):
aa:=1:bb:=1:cc:=1: #I have no idea why he allows these extra variables
A:=chrem2([bb*c,cc*q],[b,c]);
B:=chrem2([cc*a,aa*nn],[c,a]);
C:=chrem2([aa*b,bb*p],[a,b]);
mu:=gcd(A,gcd(B,C)):
A:=A/mu:B:=B/mu:C:=C/mu:
#Now get a unimodular matrix  M1  with M1*[Aa,Bb,Cc]=[1,*,*]
As :=euc3(A*a,B*b,C*c);
A1s:=l279(A*a,B*b,C*c):
#  Let M1:=matrix(3,3,[op(As),op(A1s)]):
#  Then   g:=multiply(multiply(M1,Q),transpose(M1)):
#  is equiv to our starting f, but on a 2D subspace is divisible by d=det
M2:=matrix(3,3,[op(d*As),op(A1s)]):
g1:=multiply(multiply(M2,Q),transpose(M2)):
#  g1, = his g', is then contained in f, det(g')=d^3, but is now divis by d.
#  so we let  g'''=(g1)*(1/d) ; this one is now unimodular. 
#  When we reduce it, we know what the answer will be!
M3:=ternmin(g1[1,1]/d,g1[2,2]/d,g1[3,3]/d,g1[2,3]/d,g1[3,1]/d,g1[1,2]/d):
multiply(M3,M2);
end:

#So Gauss's non-factoring alternative to  mysolve(a,b,c,n,p,q)  is  
#   gsolve(a,b,c,n,p,q):[["[2,1],"[2,2],"[2,3]],["[3,1],"[3,2],"[3,3]]];
#
#Works well on his own problem (art.294):
#gsolve(7,-15,23,3,7,6);  #yields
#[[8, -29, -23], [-105, 165, 120]]  , the last reducible to [7,11,8]
#
#But try it on, say,
#[601, 307, -101, 323, 284, 90];
#There is a huge initial swell at lemma279; then  minn  proceeds VERY slowly.
#
#A larger example  [1021, 503, -211, 602, 388, 125]  is even worse, and
#Cremona's example is essentially impossible:
#aa:=1:bb:=113922743:cc:=-310146482690273725409:
#nn:=0:p:=64039559:q:=29074923206878255640:
#gsolve(aa,bb,cc,nn,p,q);


#No witnesses handy? Try
with(numtheory):
#
hsolve:=proc(a,b,c)
gsolve(a,b,c,msqrt(-b*c,abs(a)),msqrt(-a*c,abs(b)),msqrt(-b*a,abs(c)));
[["[2,1],"[2,2],"[2,3]],["[3,1],"[3,2],"[3,3]]]:
end: