#Maple V (Release 4) code to solve  ax^2+by^2+cz^2=0  efficiently.
#From maple type
#	read `legendre.maple`:
#	mysolve( a,b,c,n,p,q );
#
#where n,p,q are 'certificates of solvability' -- see text.

#Sample problem (from Cremona preprint):
#[-113922743, 310146482690273725409, -1, 64039559, 88566846089432467791, 0]
#We get [11931641701, 7244, -7523107023591]  (with method 2), not his
#       [30106379962113, 7913, 12747947692] 
#
#Old sample: solve x^2-123452853y^2=1502226171486750244 rationally.

#Note: this is far from optimal code. In particular 
#1) We shouldn't throw away computed data when in lem4 we find a square 
#   factor of abc (I did this better in   holz()  ).
#2) Some cases of possible factors seem to occur never or only in pairs
#   (see attached data), so may indicate cases where we could tighten code;
#3) We should keep track of pairs known to be relatively prime in lem3
#4) Rather than give the simplest answers in lem1 (part2), we might give those
#   consistent with the lattice:

#Warning: the accompanying paper has undergone some reorganization and changes
#in notation, and so notation is no longer in sync with the program.
#
#Some examples and benchmarks appended at the end.
#
#Hand-crafted from recycled bits by rusin@math.niu.edu     
#Last modified 1998/08/31

with(numtheory): #for   GIgcd  .
readlib(lattice): #For  lattice  .

#Program comes from the factory with optimal (?) settings for solving
#individual numeric problems:
#method:=2:teak_k1:=true:mordell:=true:signs:=-1:verbose:=0:

#method:=0: #Use Maple's  isolve -- lame!
#method:=1: #Use simple descent only
method:=2: #Use descent + 2-dimensional lattice reduction
#method:=3: #Use 3-dimensional lattice reduction with LLL
#method:=4: #Use infallible 3-dimensional lattice reducer -- not yet available

tweak_k1:=true:  #Allow sign change in k_1 to increase cancellation
mordell:=true:   #Tweak intermediate results to stay within Holzer bound
signs:=-1:       #Test the effect of having wrong signs. (signs:=+1:disaster!)
verbose:=0:      #0 to reduce the yapping in lem4, etc.

#some global variables for experiments:
test_mode:=false:  #Allow background tests of random examples.
hitcount:=0:tabls:=array(1..30):for i to 30 do tabls[i]:=[]:od:

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:

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,expand(v1-k*v2),Q):
RETURN([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.
#with(numtheory):
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)]):
else print(`Improper call to lem1`,a,b,c,n,p,q): 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:
if (a mod u^2)<>0 then print(`lem2a error!`,a,u):RETURN():
 elif abs(u)<2 then print(`bogus lem2a`,a,b,c,u):RETURN():
 else if verbose=1 then lprint(`    Found square divisor`,u,`of a=`,a): fi:
fi:
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:
if (b mod u^2)<>0 then print(`lem2b error!`,b,u):RETURN():
 elif abs(u)<2 then print(`bogus lem2b`,a,b,c,u):RETURN():
 else if verbose=1 then lprint(`    Found square divisor`,u,`of b=`,b): fi:
fi:
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:
if (c mod u^2)<>0 then print(`lem2c error!`,c,u):RETURN():
 elif abs(u)<2 then print(`bogus lem2c`,a,b,c,u):RETURN():
 else if verbose=1 then lprint(`    Found square divisor`,u,`of c=`,c): fi:
fi:
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:
#program[A B C alpha A1 A2 A3] = text[b', a', c', gamma, n2, n3, n1]
B:=isqrt(abs(b)):
if B^2=abs(b) then 
 if B>1 then aq(a,b,c,n,p,q,1):RETURN(lem2b(a,b,c,n,p,q,B)) else
  C:=isqrt(abs(c)):
  if C^2=abs(c) then 
   if C>1 then aq(a,b,c,n,p,q,2):RETURN(lem2c(a,b,c,n,p,q,C)):
#Actually this will not occur if we begin with |B|>=|C| and |B|>1 !
   else aq(a,b,c,n,p,q,3):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:
if method=2 then 
	sol:=minv([1,k],[0,a],[abs(b),abs(c)]):
	k1:=sol[1]:k2:=sol[1]*k+sol[2]*a:
fi:
#
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 aq(a,b,c,n,p,q,4):RETURN(lem2b(a,b,c,n,p,q,betab)) fi:
betac:=gcd(beta,c):if betac>1 then aq(a,b,c,n,p,q,5):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 aq(a,b,c,n,p,q,6):RETURN(lem2c(a,b,c,n,p,q,u)):fi:
e:=gcd(d2,k2):u:=abs(d2/e): if u>1 then aq(a,b,c,n,p,q,7):RETURN(lem2b(a,b,c,n,p,q,u)):fi:
d:=gcd(b,alpha):
 e:=gcd(d,k2):if e<d then aq(a,b,c,n,p,q,8):RETURN(lem2b(a,b,c,n,p,q,d/e)):fi:
 if d>1 then aq(a,b,c,n,p,q,9):RETURN(lem2b(a,b,c,n,p,q,d)):fi:
d:=gcd(c,alpha):
 e:=gcd(d,k1):if e<d then aq(a,b,c,n,p,q,10):RETURN(lem2c(a,b,c,n,p,q,d/e)):fi:
 if d>1 then aq(a,b,c,n,p,q,11):RETURN(lem2c(a,b,c,n,p,q,d)):fi:
#Special tests if k1>1:
if gcd(k1,A)>1 then
 if verbose=1 then lprint(`    Special K...`):fi:
 sol:=lem3(a,A):if sol[4]>1 then aq(a,b,c,n,p,q,12):RETURN(lem2a(a,b,c,n,p,q,sol[4])):fi:
 if sol[5]>1 then aq(a,b,c,n,p,q,17):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 aq(a,b,c,n,p,q,18):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 aq(a,b,c,n,p,q,19):A1:=A1/e^2:alpha:=alpha*e:fi:
  e:=gcd(d,A3):if e>1 then aq(a,b,c,n,p,q,16):RETURN(lem2a(a,b,c,n,p,q,e)):fi:
 fi:
 if verbose=1 then lprint(`        ...but no square factors of abc found`):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)), 
	signs*p*k2*((a*alpha)&^(-1) mod (b/d2))], [c/d1,b/d2]);
q1:=chrem2([signs*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):
#option: try making x,y,z divisible by most of A1:
if tweak_k1=true then
 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:
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)]:
if sol<> [1,1,1] then
 aq(a,b,c,n,p,q,21):
 if verbose=1 then lprint(`    Found factor (lattice-deviation):`,sol):fi:fi:
ans:=[x,y,z]:
end:

LLLstyle:=proc(a,b,c,n,p,q)
#Use lattice reduction to find a minimal vector (=solution, by Holzer)
#Note: because of the added 2-adic condition, minimal vectors in THIS
#lattice need not be within Holzer region: example=
#[556359993597318346086961, 249440831, -160223414171, 512245395339917845762742,    52288669, 49132111880]
#Note also that this routine is sensitive to order of presentation (a,b,c)!
local x,t,rt:
if modp(a*b*c,2)=1 then
        if modp(a+b,4)=0 then x:=ODDABC(a,b,c,n,p,q):
        elif modp(b+c,4)=0 then x:=ODDABC(b,c,a,p,q,n):x:=[x[3],x[1],x[2]]:
        else x:=ODDABC(c,a,b,q,n,p):x:=[x[2],x[3],x[1]]:
        fi:
else #abc is even
        if modp(c,2)=0 then x:=EVENABC(a,b,c,n,p,q):
        elif modp(b,2)=0 then x:=EVENABC(c,a,b,q,n,p):x:=[x[2],x[3],x[1]]:
        else x:=EVENABC(b,c,a,p,q,n):x:=[x[3],x[1],x[2]]:
        fi:
fi:
if modp(x[1],2)=0 and modp(x[2],2)=0 and modp(x[3],2)=0 then 
        x:=[x[1]/2,x[2]/2,x[3]/2]:fi:
#Warning: need not be primitive! Try mysolve(1,-3221,-1216435,0,1333,  227846);
#mysolve(83, 51018545, -103, 5, 32867062, 79);
if gcd(x[1],gcd(x[2],x[3]))>1 then aq(a,b,c,n,p,q,26):fi:
if denom(x[1])>1 or denom(x[2])>1 or denom(x[3])>1 then aq(a,b,c,n,p,q,27):fi:
t:=(a*x[1]^2+b*x[2]^2+c*x[3]^2);
if t<>0 then x:=[0,0,t]:fi:
RETURN(x):
end:

ODDABC:=proc(a,b,c,n,p,q) 
local g,h,i,j,v1,v2,v3;
# Create lattice of some triples with  a u^2 + b v^2 + c w^2 = 0 mod 4abc, 
#    then find one of minimal length
# Assume  abc  is odd but a+b=0 mod 4
g:=q:g:=g/b mod c:
h:=p:h:=h/a mod b:
i:=n:i:=i/c mod a:
if i=0 then i:=1:fi: #happens if a=1
#want the lattice with  g u = v mod c, h w = u mod b, i v = w mod a
#and also  u = v mod 2, w=0 mod 2.
j:=1/(c*i) mod a:if modp(g,2)=0 then g:=g-c:fi:if modp(j,2)=1 then j:=j-a:fi:
v1:=[2*h,2*h*g*(1-c*i*j)+2*c*j,2]:v2:=[b,b*g*(1-c*i*j),0]:v3:=[0,2*a*c,0]:
v1:=SCALELATTICE(v1,v2,v3,a,b,c):
v2:="[2]:v3:=""[3]:v1:=v1[1]:
g:=a*v1[1]^2+b*v1[2]^2+c*v1[3]^2:if g=0 then RETURN(v1):fi:#it's k*4abc,|k|<3
lprint(`   LLL only hit`,g/4/a/b/c,`(not zero) with`,a,b,c,`; trying other v`);
#example: mysolve(1,-3221,-1216435,0,1333,  227846);
#Not a real example: second vector v2 is shorter than first v1!
#Also note a permutation gives a different answer! mysolve(1216435,3221,-1,227846,1333,0); has no problems like "v1 < v2".
if a*v2[1]^2+b*v2[2]^2+c*v2[3]^2=0 then aq(a,b,c,n,p,q,28): RETURN(v2):fi:
aq(a,b,c,n,p,q,29):
if a*(v2[1]+v1[1])^2+b*(v1[2]+v2[2])^2+c*(v1[3]+v2[3])^2=0 then 
        RETURN([v1[1]+v2[1],v1[2]+v2[2],v1[3]+v2[3]]):fi:
if a*(v2[1]-v1[1])^2+b*(v1[2]-v2[2])^2+c*(v1[3]-v2[3])^2=0 then 
        RETURN([v1[1]-v2[1],v1[2]-v2[2],v1[3]-v2[3]]):fi:
if a*v3[1]^2+b*v3[2]^2+c*v3[3]^2=0 then RETURN(v3):fi:
if g/4/a/b/c=-1 then 
aq(a,b,c,n,p,q,30):
#This did occur! using
#a=10^100+1243, b=10^95+151, c=-11. But only because of a loss of accuracy
#in the floating computations -- taking about Digits:=1000: seemed to work
#OK (albeit slowly).
        print(`Still failed; here's a (non-small!) solution:`);
        #from Borevich+Shafarevich p. 64
        g:=v1[1]*v1[3]+2*b*v1[2]:
        h:=v1[2]*v1[3]-2*a*v1[1]:
        i:=v1[3]^2+4*a*b:
        RETURN([g,h,i]):fi:
print(`A small solution can be found using work of B.Vallee (not coded,sorry)`);
RETURN([0,0,g/4/a/b/c]);
end:

EVENABC:=proc(a,b,c,n,p,q) 
local g,h,i,j,v1,v2,v3,B;
# Create lattice of some triples with  a u^2 + b v^2 + c w^2 = 0 mod 4abc, 
#    then find one of minimal length
# Assume  2|c but a,b are odd. Also really need c/2 to be odd:
if (c mod 4)=0 then RETURN(lem2c(a,b,c,n,p,q,2)):fi:
g:=q:g:=g/b mod c:
h:=p:h:=h/a mod b:
i:=n:i:=i/c mod a:
if i=0 then i:=1:fi: #happens if a=1
#want the lattice with  g u = v mod c, h w = u mod b, i v = w mod a
#and also  u = v mod 4, w=u*B mod 2 where ax^2+by^2+cB^2=0 mod 8
j:=1/(c*i) mod a:if modp(g,4)=3 then g:=g-c:fi:if modp(j,2)=1 then j:=j-a:fi:
  if modp(h,2)=0 then h:=h-b:fi:
v1:=[2*h,2*h*g*(1-c*i*j)+2*c*j,2]:v2:=[b,b*g*(1-c*i*j),0]:v3:=[0,2*a*c,0]:
B:=(a+b)/2 mod 2:if B=1 then v2:=[2*b+h,2*v2[2]+v1[2]/2,1]:fi:  
v1:=SCALELATTICE(v1,v2,v3,a,b,c):
v2:="[2]:v3:=""[3]:v1:=v1[1]:
g:=a*v1[1]^2+b*v1[2]^2+c*v1[3]^2:if g=0 then RETURN(v1):fi:#it's k*4abc,|k|<3
lprint(`   LLL only hit`,g/4/a/b/c,`(not zero) with`,a,b,c,`; trying other v`);
#example: EVENABC(33, -77483, 2, 8, 8360, 1);
#Not a real example: second vector v2 is shorter than first v1!
#This time, permuting abc  has no effect. 
#EVENABC(-5276187762844563,14299777,2, 454524058729886,7281560,1): same gig
#but v1 is 39% longer in squared-norm than v2 !. Also consider
#mysolve(7,2,-1,3,1,0); # versus mysolve(-1,2,7,0,1,3);
#First has v2<v1, second is fine; both give same solution {1,1,3}
if a*v2[1]^2+b*v2[2]^2+c*v2[3]^2=0 then aq(a,b,c,n,p,q,28):RETURN(v2):fi:
aq(a,b,c,n,p,q,29):
if a*(v2[1]+v1[1])^2+b*(v1[2]+v2[2])^2+c*(v1[3]+v2[3])^2=0 then 
        RETURN([v1[1]+v2[1],v1[2]+v2[2],v1[3]+v2[3]]):fi:
if a*(v2[1]-v1[1])^2+b*(v1[2]-v2[2])^2+c*(v1[3]-v2[3])^2=0 then 
        RETURN([v1[1]-v2[1],v1[2]-v2[2],v1[3]-v2[3]]):fi:
if a*v3[1]^2+b*v3[2]^2+c*v3[3]^2=0 then RETURN(v3):fi:
if g/4/a/b/c=-1 then 
aq(a,b,c,n,p,q,30):
        print(`Still failed; here's a (non-small!) solution:`);
        #from Borevich+Shafarevich p. 64
        g:=v1[1]*v1[3]+2*b*v1[2]:
        h:=v1[2]*v1[3]-2*a*v1[1]:
        i:=v1[3]^2+4*a*b:
        RETURN([g,h,i]):fi:
print(`A small solution can be found using work of B.Vallee (not coded,sorry)`);
RETURN([0,0,g/a/b/c/4]);
end:

SCALELATTICE:=proc(v1,v2,v3,a,b,c)
local e,f,g,w1,w2,w3,la:
# Find nearly-minimal vector in a 3D integer lattice
#	Note: this is pretty slow for big numbers -- probably the fault of
#	using floats; there is an integral version in Cohen's book.
#	Or maybe Maple's just lame again :-).
Digits:=5*(length(a)+length(b)+length(c)+5):
#This is a hack. Taking 3*(...) is insufficient, as the (100-95-1)-digit
#example revealed. 
e := evalf(sqrt(abs(a)));
f := evalf(sqrt(abs(b)));
g := evalf(sqrt(abs(c)));
w1:=[v1[1]*e,v1[2]*f,v1[3]*g]:
w2:=[v2[1]*e,v2[2]*f,v2[3]*g]:
w3:=[v3[1]*e,v3[2]*f,v3[3]*g]:
lprint(`Entering Maple's  lattice  command.`):
la:=lattice([w1,w2,w3]):
w1:=[round(la[1][1]/e),round(la[1][2]/f),round(la[1][3]/g)];
w2:=[round(la[2][1]/e),round(la[2][2]/f),round(la[2][3]/g)];
w3:=[round(la[3][1]/e),round(la[3][2]/f),round(la[3][3]/g)];
RETURN([w1,w2,w3]);
end:

mysolve:=proc(a,b,c,n,p,q)
#Main routine: check input, organize according to options, solve, check output
local inp,outp,ord,i,x,y,z,w,ans:
if verbose=1 then lprint(`Solving`,a,b,c,n,p,q):fi:
inp:=[a,b,c,n,p,q]:
if checkin(op(inp))=0 then RETURN([0,0,0]):fi:
if method=0 then
  ans:=eval(subs(subs({_N1=0,_N2=1,_N3=1},isolve(a*x^2+b*y^2+c*z^2)),[x,y,z])):
  fi:
if method=1 or method=2 then 
  ord:=sort.method(a,b,c):
  if (method=1 and abs(inp[ord[1]])=1) or (method=2 and 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]]: 
fi:
if method>2 then ans:=LLLstyle(op(inp)):fi:
RETURN(checkout(a,b,c,n,p,q,op(ans))):
end:

sort1:=proc(a,b,c)
#Sort so that: a1>b1>0>c1
local a1,b1,c1:
if a*b*c>0 then a1:=-a:b1:=-b:c1:=-c else a1:=a:b1:=b:c1:=c fi:
if a1<0 then
	if b1>c1 then RETURN([2,3,1]) else RETURN([3,2,1]) fi:
elif b1<0 then
	if a1>c1 then RETURN([1,3,2]) else RETURN([3,1,2]) fi:
else # c1<0 
	if b1>a1 then RETURN([2,1,3]) else RETURN([1,2,3]) fi:
fi:
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:

checkin:=proc(a,b,c,n,p,q)
#Ensure that input meets specified targets
local i:
if gcd(a,b)>1 or gcd(b,c)>1 or gcd(c,a)>1 then 
   print(`Error: need relatively prime coefficients`):RETURN(0):fi:
if a*b*c=0 then 
   print(`Error: supposed to use nonzero coefficients`):RETURN(0):fi:
if a*b>0 and a*c>0 and b*c>0 then 
   print(`Error: supposed to have different-signed coefficients`):RETURN(0):fi:
if (n^2+b*c) mod a > 0 then print(`Error: bad certificate for a`):RETURN(0):fi:
if (p^2+c*a) mod b > 0 then print(`Error: bad certificate for b`):RETURN(0):fi:
if (q^2+a*b) mod c > 0 then print(`Error: bad certificate for c`):RETURN(0):fi:
RETURN(1):
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!!`):aq(a,b,c,n,p,q,25):RETURN([0,0,0]):fi:
d:=gcd(gcd(x,y),z);if d>1 then
	if verbose=1 then lprint(`    Imprimitive solution; killing`,d):fi:
	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 
    if verbose=1 then lprint(`Normalized:`,u,`Holzer-OK`):fi:
 else aq(a,b,c,n,p,q,20):
    if verbose=1 then lprint(`Normalized:`,u,`improving...`):fi:
	if mordell=true then ans:=holz(a,b,c,n,p,q,op(ans)):fi:
 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):RETURN([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):if d>1 then aq(a,b,c,n,p,q,22) fi:
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:
 if e>1 then aq(a,b,c,n,p,q,23):fi:
 if f>1 then aq(a,b,c,n,p,q,24):fi:
 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,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:
imp_loop(a1,b1,c1,op(subs({x=n,y=p,z=q},ans)),op(subs({x=x0,y=y0,z=z0},ans))):
RETURN(subs({seq(ans[i]="[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:

##############################################################################
#   Some routines follow which are useful for testing etc.                   #
##############################################################################

trim:=proc(u) [u[1],u[2],u[3],u[4] mod u[1],u[5] mod u[2],u[6] mod u[3]]:end:

rebuild:=proc(a,b,c,x,y,z)
#Figure out an input to  mysolve  which should give output  [x,y,z]
[a,b,c,b*y/z mod a,c*z/x mod b, a*x/y mod c];
#no checking done for existence/uniqueness of this answer...
end:

#Use this to sort lists of inputs:   newlist:=sort(list,as):
as:=proc(a,b) if length(a)<length(b) then true else false fi end:

big:=proc(x,y,z)
#Compute consistent prime inputs for  mysolve  with a~x, b~y, c~z.
local a,b,c,n,p,q,ok,st:
st:=time():
c:=nextprime(z):c:=-c:
b:=nextprime(y):if (b=-c) then b:=nextprime(b):fi:
a:=x:ok:=false:
while not ok do
  a:=nextprime(a):
  lprint(`Trying prime`,a);
  ok:=true:
  if L(-a*b,-c)=-1 or L(-b*c,a)=-1 or L(-c*a,b)=-1 then ok:=false:fi:
  if a=b or a=-c then ok:=false fi:
od:
n:=msqrt(-b*c, a):
p:=msqrt(-c*a, b):
q:=msqrt(-a*b,-c):
lprint(`Setting up a problem took`,time()-st,`seconds`);
[a,b,c,n,p,q]:
end:

race:=proc(x,y,z)
#Compare various methods of solving Legendre equation
global method,tweak_k1,mordell:
local bb,st,a0,t0,b0,a1,t1,b1,a2,t2,b2,a3,t3,b3,a4,t4,b4,a1a,t1a,b1a,a1b,t1b,b1b,a1c,t1c,b1c,a1d,t1d,b1d,a2a,t2a,b2a,a2b,t2b,b2b,a2c,t2c,b2c,a2d,t2d,b2d:
bb:=big(x,y,z):
if leng([x,y,z])<50 then mordell:=false:
 method:=0:st:=time():a0:=mysolve(op(bb)):t0:=time()-st:b0:=leng(a0):fi:
method:=1:
if leng([x,y,z])<100 then
 tweak_k1:=false: mordell:=false:  
 st:=time():a1a:=mysolve(op(bb)):t1a:=time()-st:b1a:=leng(a1a):
tweak_k1:=false: mordell:=true:
 st:=time():a1b:=mysolve(op(bb)):t1b:=time()-st:b1b:=leng(a1b):
tweak_k1:=true: mordell:=false:  
 st:=time():a1c:=mysolve(op(bb)):t1c:=time()-st:b1c:=leng(a1c):
tweak_k1:=true: mordell:=true:  
 st:=time():a1d:=mysolve(op(bb)):t1d:=time()-st:b1d:=leng(a1d):fi:
method:=2:
tweak_k1:=false: mordell:=false:  
 st:=time():a2a:=mysolve(op(bb)):t2a:=time()-st:b2a:=leng(a2a):
tweak_k1:=false: mordell:=true:
 st:=time():a2b:=mysolve(op(bb)):t2b:=time()-st:b2b:=leng(a2b):
tweak_k1:=true: mordell:=false:  
 st:=time():a2c:=mysolve(op(bb)):t2c:=time()-st:b2c:=leng(a2c):
tweak_k1:=true: mordell:=true:  
 st:=time():a2d:=mysolve(op(bb)):t2d:=time()-st:b2d:=leng(a2d):
method:=3:st:=time():a3:=mysolve(op(bb)):t3:=time()-st:b3:=leng(a3):
#method:=4:st:=time():a4:=mysolve(op(bb)):t4:=time()-st:b4:=leng(a4):
lprint(`Method=0 : Time=`,t0,`; length=`,b0);
lprint(`Method=1a: Time=`,t1a,`; length=`,b1a);
lprint(`Method=1b: Time=`,t1b,`; length=`,b1b);
lprint(`Method=1c: Time=`,t1c,`; length=`,b1c);
lprint(`Method=1d: Time=`,t1d,`; length=`,b1d);
lprint(`Method=2a: Time=`,t2a,`; length=`,b2a);
lprint(`Method=2b: Time=`,t2b,`; length=`,b2b);
lprint(`Method=2c: Time=`,t2c,`; length=`,b2c);
lprint(`Method=2d: Time=`,t2d,`; length=`,b2d);
lprint(`Method=3 : Time=`,t3,`; length=`,b3);
#lprint(`Method=4: Time=`,t4,`; length=`,b4);
end:

leng:=proc(u) local v: 
v:=u[1]*u[2]*u[3]: if v=0 then 1 else evalf(log(abs(v))/log(10)) fi end:

aq:=proc(a,b,c,n,p,q,zz)
#Tally interesting occurences within algorithms
local i:
global hitcount,tabls:
tabls[zz]:=[op(tabls[zz]),[a,b,c,n,p,q]]:
hitcount:=hitcount+1:
if test_mode=true then
  if hitcount>100 then 
	save(tabls,TABLS):hitcount:=0:fi:fi:
end:

do_me:=proc(r,s,t)     
#Cause a lot of examples to be run up to [a,b,c]=[10^12r, 10^12s, 10^12t]
local bb,aa,cc,outp,bigt,solvt,st:
cc:=0:outp:=[]:bigt:=0:solvt:=0:
while 1=1 do
	st:=time():
	bb:=big(rand()^r,rand()^s,rand()^t):
	bigt:=bigt+time()-st:
	st:=time():
	aa:=mysolve(op(bb)):
	solvt:=solvt+time()-st:
	#outp:=[op(outp),[op(aa),op(bb)]]:
	outp:=[op(outp),[bigt,solvt]]:
	cc:=cc+1:
	if cc mod 100 = 0 then save(outp,OUTP):fi:
od:
end:

#Next line allows command-line processing: from UNIX prompt just type
#      maple  < legendre.maple > /dev/null &
#and many random examples will be generated and studied as above
#
 if test_mode=true then do_me(3,3,1); fi:

##############################################################################
fi;  # that's an intentional error to stop the reading of this file.
##############################################################################

Interesting samples of (a,b,c,n,p,q) caught by the various checks (method:=2).

1: [7, -4, 1, 5, -3, 0]
   [-5, 4, 1, 4, 3, 0]
   [61, -4, 3, 16, 30, -2]
2: none possible with method=2 if abc sorted by magnitude
3: [8, 1, -1, 5, 0, 0]
   [-5, 1, 1, -3, 0, 0]
   [525973, 1, -1, 150277, 0, 0]
   (etc. ALL c=+-1)
4: [22729, 32, -1, -6722, 3, 0]
   [-4097, 27, 5, 3554, 17, 2]
   [21127, 48, -7, 13890, 47, 6]
5: [8779, 17, -12, 496, -4, -4358749]
   [15182006953, -40426, 9, 10023777100, 14805, 1]
6: [-4103, 23, 8, 2943, 16, 7]
   [-3140477, 43, 25, 2136947, 17, 19]
   [-239743087, 78, 25, 238041788, 1, 6]
7: [-134593, 56, 1, 60372, -9, 0]
   [-47287, 8, 7, 27215, 1, 4]
   [-4431282451, 4719, 7, 582020523, 710, 6]
8: [-336649, 300, 1, 320403, 157, 0]
   [51018545, -412, 83, 14715579, 79, 10]
   [134873197, 1791, -13, 124380370, 838, 8]
9: none found
10: [9400805511, 2437, -4, 48738319, 1068, 1]
    [2607882701, -17205, 4, 2216847176, 4106, 3]
11: none found
12: [153, -13, 1, 25, 9, 0]
    [-15152, 5, 3, 12535, 1, 1]
    [4293, 13, -10, 884, 11, 9]
    [32432, -15, 7, 29773, 1, 1]
13: unused
14: unused
15: unused
16: none found
17: [7185786, -1465, 1, 5630881, 268, 0]
    [21328386, 4537, -1, 20187347, 1468, 0]
    [-30940374, 1871, 1, 25555033, 384, 0]
    [29239123130, -165321, 1, 6620397451, 5579, 0]
    [86807326430, 57649, -1, 21462227647, 9126, 0]
    [-17194283526, 62471, 1, 14882593087, 21917, 0]
    [130017558910, 31761, -1, 119033273479, 6361, 0]
    [188060224470, -31489, 1, 68317415963, 14697, 0]
    [-2568561397757146, 1015687, 9, 61755002965489, 70774, 8]
18: Note: all examples found imply additional examples of 17, too (above are
	17-ers NOT also flagging #18). Some samples:
    [6819, 29, -5, 991, 7, 3]
    [9987, -91, 1, 5753, 32, 0]
    [-9476, 319, 1, 8563, 246, 0]
    [-242324, 61, 43, 56873, 41, 27]
    [132963927, -1145, 2, 101950561, 336, 1]
19: none found
20: [7,2,-1,3,1,0]
    [5, -83, -2, 582, 33, 285]
    [257, -116044, -1, 67043400, 38909, 0]
    [1027986085027861747937, 722974121783, -604305613931, 97844783073873996625, 394665269920, 282075028651]
    [21458238196506826611317, 155590763477, -429392673713, 783370648399676997236, 92269155657, 136629727139]
    [293911884200387681576093, 106205605559, -687223946819, 270630693817753504680542, 1203796437, 54028739034]
    [206207486824798778186641, 769034410087, -80703111317, 169619646075278274214379, 489267596820, 71563083678]
    (Only occasionally need to loop more than once!. But these last two are
    examples to show this can be necessary. Last even spooks LLL !)
21: [79, -8, 1, 61, 843, 0]
    [125965, -134, 1, 78828, 53, 0]
    [-163184, 223, 1, 39551, 174, 0]
    [-1115479181, 2065, 61, 998731466, 674, 24]
    [88727009405, -50601, 1, 47011375381, 18361, 0]
    [263422353799, -497480, 1, 76811143383, 115159, 0]
   [62063469603854588654579,68044230041,-1,766592503106114120003,39384669929,0]
    (many. Usually: close to abc=1,-1,c; c=+-1;factors small; only one lattice 
    congruence wrong; deducible from descended violation; etc. 
    Here I've listed some more interesting failures.)
22: [1, 7535987, -52, 0, 3826534, 49]
    [19, 4776762553, -6836, 14, 395905266, 2005]
    [1, 52753834808, -208881, 0, 25716985333, 95464]
    [221916530122928255, 1, -370095536, 108860950807351256, 0, 232494953]
23: Only (?) occur when  24  also does (but not vice versa).
    [555143, 1, -1737, 15788, 0, 1231]
    [1, 23121989, -2313, 0, 19425797, 586]
    [14299777, 2, -5276187762844563, 7281560, 1, 454524058729886]
24: [1, 316, -5, 0, 99, 3]
    [388, 152441, -1, 259, 29878, 0]
    [831, 1, -92742160, 464, 0, 46170687]
    [1, 29526647, -1956, 0, 13408522, 1885]
    [7, 4715959977, -21016, 4, 1325206462, 11053]
    [6404148857786219, 1, -3120332, 5132340781795971, 0, 143151]

Frequency (evidently 1705 cases)
> for i to 25 do lprint(i,nops(tabls[i]));od;
1   501
2   0
3   1705
4   167
5   4             10 in only 1200 big rands
6   1              2
7   64            46
8   11            20
9   0
10   1             4
11   0
12   95           98
13   0
14   0
15   0
16   0
17   54           58
18   47           50
19   0
20   958
21   1729
22   3
23   6
24   108         140
25   0

For  do_me(8,5,2):
1   335
2   0
3   1121
4   172
5   6
6   3
7   67
8   27
9   0
10   1
11   0
12   115
13   0
14   0
15   0
16   0
17   85
18   70
19   0
20   933
21   1992
22   7
23   11
24   206

Averages of 1700 trials of the form big(rand^2,rand,rand):
	1.061784118 sec to create ( big() )
	0.3761 sec to solve ( mysolve() )  [method 3 average: 6sec]
Averages of 1200 trials of the form big(rand^3,rand^2,rand):
	1.9542 sec to create ( big() )
	0.5032 sec to solve ( mysolve() )
AAverages of 3600 trials of the form do_me(3,3,1):
	2.2436 sec
	0.7179 sec
Averages of 100 trials of the form big(rand^8,rand^5,rand^2):
	18.5349 sec to create ( big() )
	 1.1943 sec to solve ( mysolve() )
	(for 1100 trials these numbers were [23.39942818, 1.384026364])

Race results:
race(10^100+1242,10^95+150,10);#next 3 integers are compatible primes.
Method=1a: Time=   10.460   ; length=   220.0517970
Method=1b: Time=   230.900   ; length=   195.1141345
Method=1c: Time=   11.330   ; length=   214.4273152
Method=1d: Time=   11.720   ; length=   194.8751826
Method=2a: Time=   1.210   ; length=   208.8157116
Method=2b: Time=   1.170   ; length=   195.1141345
Method=2c: Time=   1.190   ; length=   196.1624390
Method=2d: Time=   1.280   ; length=   195.1141345
Method=3 : Time=   15.240   ; length=   1  [only hit  -1 ; BS trick invalid?]

race(10^150+1458, 10^100+266,10^50+150);
Method=1a: Time=   t1a   ; length=   b1a
Method=1b: Time=   1584.520   ; length=   296.9701225
Method=1c: Time=   30.750   ; length=   319.5262711
Method=1d: Time=   30.400   ; length=   296.9701225
Method=2a: Time=   3.320   ; length=   296.9701225
Method=2b: Time=   3.150   ; length=   296.9701225
Method=2c: Time=   3.190   ; length=   296.9701225
Method=2d: Time=   3.040   ; length=   296.9701225

st:=time():bb:=big(10^300,10^200,10^100);time()-st;
	time = 3003.270

bb := [10000000000000000000000000000000000000000000000000000000000000000000\
00000000000000000000000000000000000000000000000000000000000000000000000\
00000000000000000000000000000000000000000000000000000000000000000000000\
00000000000000000000000000000000000000000000000000000000000000000000000\
00000000000000009559, 10000000000000000000000000000000000000000000000000\
00000000000000000000000000000000000000000000000000000000000000000000000\
00000000000000000000000000000000000000000000000000000000000000000000000\
000000357, -100000000000000000000000000000000000000000000000000000000000\
00000000000000000000000000000000000000267, 80171399393609977980132855940\
17716334629251822906657500894399141590357757863311873036803998648122405\
63173807584536542031017197851577452186119513340576362591365283403643335\
18057262171819564245602293369015687921239062523424040311203770822669673\
9726283539921459515329307052593554238048263654497759457654, 614961432997\
43995890445687791180511148274988143934636252819628494836689651438752547\
93040897285508680771611443469717164146144461962349179520660517642948613\
459528986289871669120982647790175966826032206, 6419270008326324782117042\
60624775611960310424188820665671462777663688494030548668562567762400539\
8289]

> method;       2
> tweak_k1;     true
> mordell;      true
> st:=time():aa:=mysolve(op(bb)):tt:=time()-st;

[several calls with long reports deleted]
    Special K...
        ...but no square factors of abc found
Solving   2366841662820213134   -6454846499286591487550656043095201381   -1   
1897589059675679063   1524237756364834510416061973151271222   0
Solving   109594523   2366841662820213134   -1   
781184565206609718181508123640   1793919991421482803   0
Solving   3197   109594523   -1   311222047217544   -2887627   0
Solving   1   3197   -1   0   -3196   0
Normalized:   [.1768596177e-1, 0, .1768596177e-1]   Holzer-OK
Normalized:   [.7613141399e-1, .1768596177e-1, .7815871955e-1]   Holzer-OK
Normalized:   [.7598136535e-1, .7613141399e-1, .1075600301]   Holzer-OK
Normalized:   [.1646252768, .7598136535e-1, .1460421648]   Holzer-OK
Normalized:   [.3839823961e-1, .1646252768, .1600845306]   Holzer-OK
    Found factor (lattice-deviation):   [3, OK, OK]
Normalized:   [.4033614634e-1, .1279941320e-1, .4231819563e-1]   Holzer-OK
    Found factor (lattice-deviation):   [OK, 3, OK]
Normalized:   [.1909980238e-1, .4033614634e-1, .4462966675e-1]   Holzer-OK
Normalized:   [.5729940714e-1, .1271863484e-2, .5731352104e-1]   Holzer-OK
                                 tt := 19.880

This is amazingly fast! After nearly an hour's search for a problem, it's
solved with 20 seconds of CPU time!

method1, with tweak_k1 and mordell still true: tt := 156.490