#This file contains Maple V (Releases 2..5) code borrowed (with only minor
#adaptations) with thanks from Dave Rusin, to solve  ax^2+by^2+cz^2=0  
#efficiently, based on his collaborative work with John Cremona. See
#    http://www.math.niu.edu/~rusin/research-math/legendre/
#    http://www.maths.nott.ac.uk/personal/jec/

chrem2:=proc(a,b) local u;
#Maple idiocy: chrem([a1,a2],[b1,1]) ==> ERROR, so we substitute
  u:=`mod`(a[1],b[1]):
  u+b[1]*`mod`((a[2]-u)/b[1],b[2])
end:

minv:=proc(v1,v2,Q)
#Express the smallest linear combo (relative to Q) of v1 and v2 as x v1 + y v2.
#  This has come to be the most time-critical loop!
local L1,L2,IP,sol,k:
L1:=Q[1]*v1[1]^2+Q[2]*v1[2]^2:
L2:=Q[1]*v2[1]^2+Q[2]*v2[2]^2:
if L2>L1 then sol:=minv(v2,v1,Q):return [sol[2],sol[1]] fi:
IP:=Q[1]*v1[1]*v2[1]+Q[2]*v1[2]*v2[2]:
k:=round(IP/L2):
if k=0 then return [0,1] fi:
if k=1 and IP/L2=1/2 then return [0,1] fi:
sol:=minv(v2,[v1[1]-k*v2[1],v1[2]-k*v2[2]],Q):
[sol[2],sol[1]-k*sol[2]]
end:

lem1:=proc(a,b,c,n,p,q)
#Find a solution directly for cases with |ab|=1 or... 
local u,a1,b1,c1:
if a*b=-1 then return [1,1,0] fi:
if b*c=-1 then return [0,1,1] fi:
if c*a=-1 then return [1,0,1] fi:
#So assume ab=1 etc.
if a*b=1 then
 if a=-1 then c1:=-c else c1:=c: fi:
 u:=GIgcd(q+I,c1):
 return [Re(u),Im(u),1]
elif b*c=1 then
 if b=-1 then a1:=-a else a1:=a: fi:
 u:=GIgcd(n+I,a1):
 return [1,Re(u),Im(u)]
elif a*c=1 then
 if a=-1 then b1:=-b else b1:=b: fi:
 u:=GIgcd(n+I,b1):
 return [Im(u),1,Re(u)]
fi
end:

lem2a:=proc(a,b,c,n,p,q,u)
#Peel off a known square factor  u  from coefficient  a
local p1,q1,sol:
p1:=p*(u&^(-1) mod b) mod b;
q1:=q*(u&^(-1) mod c) mod c;
sol:=mysolve(a/u^2,b,c,n,p1,q1):
p1:=gcd(sol[1],u):
[sol[1]/p1,u*sol[2]/p1,u*sol[3]/p1]
end:

lem2b:=proc(a,b,c,n,p,q,u)
#Peel off a known square factor  u  from coefficient  b
local n1,q1,sol:
n1:=n*(u&^(-1) mod a) mod a;
q1:=q*(u&^(-1) mod c) mod c;
sol:=mysolve(a,b/u^2,c,n1,p,q1):
q1:=gcd(sol[2],u):
[u*sol[1]/q1,sol[2]/q1,u*sol[3]/q1]
end:

lem2c:=proc(a,b,c,n,p,q,u)
#Peel off a known square factor  u  from coefficient  c
local n1,p1,sol:
n1:=n*(u&^(-1) mod a) mod a;
p1:=p*(u&^(-1) mod b) mod b;
sol:=mysolve(a,b,c/u^2,n1,p1,q):
n1:=gcd(sol[3],u):
[u*sol[1]/n1,u*sol[2]/n1,sol[3]/n1]
end:

lem3:=proc(a,b)
#Decompose a, b into relatively prime parts and squares.
local m1,m2,m12,b1,b2,u,check_em:
#
check_em:=proc(x,y,z)
#See if x,y,z relatively prime else return corrections to use as  u
local d,ans:
d:= gcd(x,y):
if d>1 then ans:=[1/d,1/d,d,1,1]:
else
  d:=gcd(x,z):
  if d>1 then ans:=[1/d,d,1/d,d,1]:
  else
    d:=gcd(y,z):
    if d>1 then ans:=[d,1/d,1/d,1,d]:
    else ans:=[]:
    fi:
  fi:
fi:
ans
end:
#
m1:=a:
m2:=b:
m12:=1:
b1:=1:
b2:=1:
u:=[1,1,1,1,1]:
while u<>[] do
    m1:=m1*u[1]:
    m2:=m2*u[2]:
    m12:=m12*u[3]:
    b1:=b1*u[4]:
    b2:=b2*u[5]:
    u:=check_em(m1,m2,m12):
    od:
[m1,m2,m12,b1,b2]
end:

lem4:=proc(a,b,c,n,p,q)
#Main descent: Either
# Find a square factor of  abc  and call lem2,
# Or, find  bc=+-1 and call lem1, or
# Or, compute a smaller (A,B,C,n1,p1,q1) and call lem4 (again).
local k1,k2,t,k,A,B,C,alphA,betA,d1,d2,betAb,betAc,d,e,u,v,n1,p1,q1,sol,\
x1,y1,z1,x,y,z,ans,A1,A2,A3,x2,y2,z2:
B:=isqrt(abs(b)):
if B^2=abs(b) then 
 if B>1 then return lem2b(a,b,c,n,p,q,B) else
  C:=isqrt(abs(c)):
  if C^2=abs(c) then 
   if C>1 then return lem2c(a,b,c,n,p,q,C)
#Actually this will not occur if we begin with |B|>=|C| and |B|>1 !
   else return lem1(a,b,c,n,p,q) fi:
  fi:
 fi:
fi:
#
k:=n*(c&^(-1) mod a) mod a;
if k>abs(a)/2 then k:=k-abs(a):fi:
k1:=1:k2:=k:
sol:=minv([1,k],[0,a],[abs(b),abs(c)]):
k1:=sol[1]:k2:=sol[1]*k+sol[2]*a:
#
t:=(b*k1^2+c*k2^2)/a:
sol:=lem3(t,b*c):
A:=sol[1]:B:=sol[2]:C:=sol[3]:alphA:=sol[4]:betA:=sol[5]:
if C>0 then C:=-C:A:=-A:B:=-B: fi:
#Eight chances to find a square factor of abc:
betAb:=gcd(betA,b):if betAb>1 then return lem2b(a,b,c,n,p,q,betAb) fi:
betAc:=gcd(betA,c):if betAc>1 then return lem2c(a,b,c,n,p,q,betAc) fi:
d1:=gcd(C,c):d2:=C/d1:
e:=gcd(d1,k1):u:=abs(d1/e): if u>1 then return lem2c(a,b,c,n,p,q,u) fi:
e:=gcd(d2,k2):u:=abs(d2/e): if u>1 then return lem2b(a,b,c,n,p,q,u) fi:
d:=gcd(b,alphA):
 e:=gcd(d,k2):if e<d then return lem2b(a,b,c,n,p,q,d/e) fi:
 if d>1 then return lem2b(a,b,c,n,p,q,d) fi:
d:=gcd(c,alphA):
 e:=gcd(d,k1):if e<d then return lem2c(a,b,c,n,p,q,d/e) fi:
 if d>1 then return lem2c(a,b,c,n,p,q,d) fi:
#Special tests if k1>1:
if gcd(k1,A)>1 then
 sol:=lem3(a,A):if sol[4]>1 then return lem2a(a,b,c,n,p,q,sol[4]) fi:
 if sol[5]>1 then A:=A/sol[5]^2:alphA:=alphA*sol[5]:fi:
 A1:=sol[2]:A2:=sol[3]:A3:=sol[1]:
 d:=gcd(k1,alphA):if d>1 then k1:=k1/d:k2:=k2/d:alphA:=alphA/d:fi:
 d:=gcd(k1,k2):
 if d>1 then 
  e:=gcd(d,A1):if e>1 then A1:=A1/e^2:alphA:=alphA*e:fi:
  e:=gcd(d,A3):if e>1 then return lem2a(a,b,c,n,p,q,e) fi:
 fi:
else A1:=A:A2:=1:fi:
#Enough tests! Find the new certificates
n1:=chrem2([c*k2*(k1&^(-1) mod A1),n ],[A1,A2]):
p1:=chrem2([q*k1*((a*alphA)&^(-1) mod (c/d1)), 
	-p*k2*((a*alphA)&^(-1) mod (b/d2))], [c/d1,b/d2]);
q1:=chrem2([-p*A*alphA*(k1&^(-1) mod d2), q*A*alphA*(k2&^(-1) mod d1)],
	[d2,d1]);
#Descend & climb back:
sol:=mysolve(A,B,C,n1,p1,q1):
x1:=sol[1]:y1:=sol[2]:z1:=sol[3]:
x:=-alphA*A*x1:
y:=(c/d1)*(k2/d2)*y1+k1*z1;
z:=(b/d2)*(k1/d1)*y1-k2*z1:
d:=gcd(y,z):x1:=gcd(d,x):
#try making x,y,z divisible by most of A1:
y2:=(c/d1)*(k2/d2)*y1-k1*z1;
z2:=-(b/d2)*(k1/d1)*y1-k2*z1:
e:=gcd(z2,y2):x2:=gcd(e,x):
if x2>x1 then y:=y2:z:=z2:x1:=x2:d:=e:fi:
if x1>1 then x:=x/x1:y:=y/x1:z:=z/x1:d:=d/x1:fi:
#if d>1, then d^2 | a ...
sol:=[should(abs(a),n*z,b*y),should(abs(b),p*x,c*z),should(abs(c),q*y,a*x)]:
ans:=[x,y,z]
end:

mysolve:=proc(a,b,c)
#Main routine: 
local inp,outp,ord,i,ans,n,p,q:
if nargs=6 then n:=args[4]:p:=args[5]:q:=args[6]:
else n:=msqrt(-b*c,abs(a)):p:=msqrt(-a*c,abs(b)):q:=msqrt(-b*a,abs(c)):fi:
inp:=[a,b,c,n,p,q]:
ord:=sort2(a,b,c):
if abs(inp[ord[2]])=1
  then outp:=lem1(seq(inp[ord[i]],i=1..3),seq(inp[ord[i]+3],i=1..3)):
else outp:=lem4(seq(inp[ord[i]],i=1..3),seq(inp[ord[i]+3],i=1..3)): fi:
ans:=array(1..3):
for i to 3 do ans[ord[i]]:=outp[i]: od:
ans:=[ans[1],ans[2],ans[3]]: 
checkout(a,b,c,n,p,q,op(ans))
end:

sort2:=proc(a,b,c)
#Sort so that: |a|>|b|>|c|
local a1,b1,c1:
a1:=abs(a):b1:=abs(b):c1:=abs(c) :
if a1>b1 then
	if a1>c1 then 
		if b1>c1 then return [1,2,3] else return [1,3,2] fi:
	else return [3,1,2] fi:
else
	if a1<c1 then 
		if b1<c1 then return [3,2,1] else return [2,3,1] fi:
	else return [2,1,3] fi:
fi
end:

checkout:=proc(a,b,c,n,p,q,x,y,z)
#Check that purported solution is as desired -- if not, fix if possible
local u,d,ans:
ans:=[x,y,z]:
if a*x^2+b*y^2+c*z^2<>0 then 
	print(Bad_solution):return [0,0,0] fi:
d:=gcd(gcd(x,y),z);if d>1 then
	ans:=[x/d,y/d,z/d]:fi:
u:=evalf([sqrt(ans[1]^2/abs(b*c)),sqrt(ans[2]^2/abs(a*c)),sqrt(ans[3]^2/abs(b*a))]):
if u[1]<=1 and u[2]<=1 and u[3]<=1 then 
    ans:=holz(a,b,c,n,p,q,op(ans)):
fi:
ans
end:

euc:=proc(a,b,c)
#Euclidean algorithm: return x,y where a x - b y = c  (gcd(a,b)|c is ASSUMED)
local s,k:
if a=0 then return [0,-c/b] fi:
if b=0 then return [c/a,0] fi:
if a<0 then s:=euc(-a,b,c):return [-s[1],s[2]] fi:
if b<0 then s:=euc(a,-b,c):return [s[1],-s[2]] fi:
if a<b then s:=euc(b,a,c):return [-s[2],-s[1]] fi:
k:=round(a/b): s:=euc(a-b*k,b,c):[s[1],s[2]+k*s[1]]
end:

improve:=proc(a,b,c,x,y,z)
#Use trick found in [Mordell] to make a solution just a little better
local k,e,u,v,w,r,s,par:
if c mod 2 = 0 then
 k:=c/2:
 e:=euc(y,x,k):u:=e[1]:v:=e[2]:
 w:=round(-(a*u*x+b*v*y)/c/z):
else
 k:=c*2:
 e:=euc(y,x,c):u:=e[1]:v:=e[2]:
 par:=(a*u+b*v) mod 2:
 w:=round(-(a*u*x+b*v*y)/c/z):
 if (w-par) mod 2 = 1 then 
	if w>(-(a*u*x+b*v*y)/c/z) then w:=w-1 else w:=w+1:fi:
 fi:
fi:
r:=a*u^2+b*v^2+c*w^2:
s:=a*u*x+b*v*y+c*w*z:
[(x*r-2*u*s)/k,(y*r-2*v*s)/k,(z*r-2*w*s)/k]
end:

imp_loop:=proc(a,b,c,n,p,q,x,y,z)
#Loop over improvements in [Mordell] until Holzer bound attained.
#  Here we ASSUME gcd(x,y,z)=1  (so gcd(x,y)^2 | c).
local s,x0,y0,z0,prev,d,e,f,c0:
d:=gcd(y,x):
prev:=abs(z)+1:
x0:=x/d:y0:=y/d:z0:=z:c0:=c/d^2:
while z0^2>(a*b) do 
 if abs(z0)>=prev then 
   print(imp_loop_problem,a,b,c,x,y,z):return [x0,y0,z0]
 fi:
 s:=improve(a,b,c0,x0,y0,z0):prev:=abs(z0):
 e:=gcd(s[1],s[2]):f:=gcd(e,s[3]):e:=e/f:
 x0:=s[1]/e/f:y0:=s[2]/e/f:z0:=s[3]/f:c0:=c0/e^2:d:=d*e:
 od:
e:=gcd(d*gcd(x0,y0),z0):
s:=[x0*d/e,y0*d/e,z0/e]
end:

holz:=proc(a,b,c,n,p,q,x0,y0,z0)
#Adjust a big solution to get one in Holzer region. ASSUME gcd(x0,y0,z0)=1.
local a1,b1,c1,ans,i,w,x,y,z:
if a*b*c>0 then a1:=-a:b1:=-b:c1:=-c else a1:=a:b1:=b:c1:=c: fi:
ans:=[x,y,z]:
if a1<0 then w:=c1:c1:=a1:a1:=w: ans:=[ans[3],ans[2],ans[1]] fi:
if b1<0 then w:=c1:c1:=b1:b1:=w: ans:=[ans[1],ans[3],ans[2]] fi:
w:=imp_loop(a1,b1,c1,op(subs({x=n,y=p,z=q},ans)),op(subs({x=x0,y=y0,z=z0},ans))):
subs({seq(ans[i]=w[i],i=1..3)},[x,y,z])
end:

should:=proc(a,b,c)
#Check to see if  b is congruent to  +- c  mod a  (assumed positive!)
#  if not, how much to divide  a  by to ensure that congruence holds?
local u,v: u:=gcd(a,b-c):v:=gcd(a,b+c): a/max(u,v)
end: