#Here is a suite of Maple procedures for performing descent on elliptic 
#curves of the form  y^2=x(x^2+Ax+B). The key novel feature is the ability
#to check homogeneous spaces y^2=d1x^4+Ax^2+(B/d1) which pass the usual
#2-descent tests to see if they pass a second 2-descent test. The tests
#should run efficiently if a few initial parameters can be factored.
#
#Read this file into maple and type
#      MAXRANK(d1,A,B);
#You'll get a number which is the maximum rank possible for the curve
#consistent with the descent calculations. (In practice it equals the rank.)
#Some more comments are available by typing "HELP();"
#
#Supporting information describing the routines and why they work will
#eventually appear at my web site 
#(http://www.math.niu.edu/~rusin/research-math/descent/)
#The methods are similar to those in papers of Merriman,Siksek,Smart;
#Birch,Swinnerton-Dyer; Pollard, Schnorr; and to programs of 
#Connell, and Cremona, but seem to be relatively unknown.
#
#Dave Rusin, rusin@math.niu.edu
dat:=`1997/05/02`:

##################### INITIAL COMMANDS ############################

lprint(`Please excuse a dozen lines of gibberish...`);
with(numtheory):#These lines can be streamlined for space-tight jobs
with(padic):
with(linalg): 
with(combinat):#we use only   powerset
readlib(ifactors):
readlib(issqr):
readlib(iroot):
readlib(symmdiff):
readlib(realroot):
readlib(lattice):
isolve(x^2+y^2=2*z^2):print(`numtheory/msqrt`):print(_msqrtprime):print(_msqrtpk):
#Dont ask me why but that is the easiest way to make the proc's available!

read `prims.m`: #a file stating KNOWNPR:={2,3,[your favorites]}: 
#insert any primes worth knowing. 
if KNOWNPR='KNOWNPR' then read(prims):fi:#sometimes I use a text file
if KNOWNPR='KNOWNPR' then KNOWNPR:={2}:fi:save(KNOWNPR,`prims.m`):
UNKNOWN:={}:#unfactored numbers, fatal or not. See FACTORING section.
FATALSET:=[]:#composites to check with factorization program for future runs

yaklevel:=4:#values -1 0 1 ... 10  have meaning -- see PRT();
 #If used to create data for a larger program, then  yaklevel:=0 is appropriate
printlevel:=3:#I like to have maple be more informative with error messages...
factorproc:=`external`: #method to be used for large integers
 #other options include  easy  maple (=Morrison-Brillhart) lenstra  etc.
FACTLIM:=2^96:#What constitutes "small" for factoring choice?
isprime(1093^2):=false: #maple gets this wrong!
isprime(3511^2):=false: #it would thus mess up factoring, too
#gc(500000):# halves rate of garb. collection 
gc(0): # preserves rate, but stops displaying "bytes used" etc
FACTTIME:=0:#time spent factoring
HASSETIME:=0:#time spent collecting divisors allowing solvable conics
NORMTIME:=0:#time spent finding points on conics
LSTIME:=0:#time spent with Siksek's local solvability criteria for quartics
QTIME:=0:#time spent minimizing and reducing quartics
 #these time counts overlap
EXPERIMENTAL:=false:#no experiments currently in progress anyway.
WANTSEARCH:=true:#save input for  ratpoint?  Note: overwrites old files!
JUSTONE:=false:#"true" means find one descendant quartic for each d1

#global variables the user is not supposed to think about:
SNAG:=0: LIKELIES:=[2]:BADSET:={}:AMBIG:=0:FATAL:=false:TORSION:=1:SEARCHES:=0:

lprint(`Carry out 1st and 2nd 2-descent to bound the rank of the curve y^2=x(x^2+Ax+B),`);
lprint(`    by typing    MAXRANK(A,B);   Type   HELP();  for more info.`);
lprint(`Reading in the procedures will take just a minute.`);
lprint();

#######################  FRONT END 'N' STUFF  ##############################
# These are the only routines the user should call; the others below often #
# _assume_ that type-checking is done, that global variables are set,      #
# that GCD's=1, etc.                                                       #
############################################################################

#Different entry points according to how broad the user's problem is. These
#are all standard reductions. All lead to the guts of the program (the
#higher descents) through the routine  GUTS.

MAXRANK:=proc(A,d1d2) 
local ti,AA,dd,f,i,W,native,maxim,found;
global EXPERIMENTAL:
#
#Estimate the ranks both of  E  and the 2-isogenous curve  E'.
#
if nargs=0 then print(`Type  HELP();  for help.`);RETURN():fi:
LRT(2,`We compute an upper bound on the rank of the curve`);
PRT(2,Y^2=X*(X^2+A*X+d1d2));
#Secret Code (TM) to try the latest experimental routines (Shh!):
if nargs>2 then 
	EXPERIMENTAL:=true:WARN(`Using experimental procedures.`):
else EXPERIMENTAL:=false:fi:
ti:=time();
PRETTIFY(A,d1d2):AA:="[1]:dd:=""[2]:
if TORSION=0 then RETURN(0):fi:#give up  -- curve is singular.
LRT(2,`We first look for generators in E/phi(E').`);
native:=ONECURV(AA,dd):
LRT(2,` _,-'"``-._,-'"``-._,-'"``-._,-'"``-._,-'"``-._,-'"``-._,-'"``-._,-'"``-._`);
LRT(2,`We now look for generators in phi(E')/(2E).`);
PRETTIFY(-2*AA,AA^2-4*dd):#just to set TORSION I guess
maxim:=ONECURV(-2*AA,AA^2-4*dd):
LRT(2,` _,-'"``-._,-'"``-._,-'"``-._,-'"``-._,-'"``-._,-'"``-._,-'"``-._,-'"``-._`);
LRT(2,` `):
LRT(2,`Rank of the original curve is at most`,native+maxim);
LRT(3,`This curve required`,time()-ti,`seconds to analyze.`);
LRT(3,` `);
RETURN(native+maxim):
end:

ONECURV:=proc(A,d1d2) #extra args=x-coords of any known generators.
local goods,AA,dd,f,P,i,W,new,W0,found,found2,s,minim,maxim;
#
#Compute the free rank of  E/im(phi). 
# Note: we should use group structure to: solve Norm eqns, determine loc. solv,
# etc.; we could change basis whenever we get a FATAL and thus keep going.
# Not yet implemented.
#
if nargs=0 then
	lprint(`ONECURV(A,B): Compute the number of free generators of E:y^2=x(x^2+Ax+B) `);
	lprint(`which are not in the image of the obvious 2-isogeny. Any additional`);
	lprint(`arguments assumed to be x-coordinates of points on  E (speeds things up)`);
	lprint(`For more info try  HELP(): `); RETURN():
fi:
LRT(2,`We bound the free rank of  E/im(phi)  where  E  is the curve`):
PRT(2,Y^2=X*(X^2+A*X+d1d2));
PRETTIFY(A,d1d2):AA:="[1]:dd:=""[2]:
if TORSION=0 then RETURN(0);fi:
goods:={}: minim:=0:
W0:=ODDINT(dd):P:=max(op(W0)):#use the 2-torsion element to toss largest prime
if TORSION>1 then 
 f:=ODDINT((-AA+sqrt(AA^2-4*dd))*2):s:=1:for i in f do s:=s*i:od:
 if member(P,f) then f:=symmdiff(f,W0):fi:
 found:={1,s}:found2:={{},f}:
else found:={1}: found2:={{}}: 
fi:
if (nargs>2) or (TORSION>1) then 
 #allow user to enter x-coordinates of known generators to save time!
 if nargs>2 then LRT(4,`You have entered`,nargs-2,`points already known.`):fi:
 for i from 3 to nargs do 
   f:=args[i]:
   if not type(sqrt(f*(f^2+A*f+d1d2)),rational) then 
	ERR(f,`isn't an x-coordinate of a point.`):
   else goods:={op(goods),f}: fi:
 od:#Now each of these good points doubles the number of h-sp with points:
 for i in goods do 
  W:=ODDINT(i); if member(P,W) then W:=symmdiff(W,W0):fi:
  new:={}:for s in found2 do new:={op(new),s,symmdiff(s,W)}:od:
  found2:=new:# found2 just doubled in size.
  od:
 i:=0:while nops(found2)>2^i do i:=i+1:od:minim:=i:
 for f in found2 do s:=1:for i in f do s:=s*i:od: found:={op(found), s}:od:
fi:
#OK, now let's list the d1's (i.e. the homog.spp.) which can have points:
W:=FIND_D1S(AA,dd) minus found: #really should be a _quotient_.
found:=found union MANYD1(W,AA,dd) union {1}; #that's the crux of the procedure
#find the largest sub_group_ of  "found"....
i:=0:while nops(found)>2^i do i:=i+1:od:maxim:=i:
if TORSION=2 then maxim:=maxim-1:minim:=minim-1:fi:
 #our counts have modded out by the special element of order 2 only.
if maxim=0 then LRT(2,`The curve has rank zero.`):RETURN(0):fi:
if minim=maxim then LRT(2,`All`,minim,` generators have been found (up to odd index)`);fi:
if nargs>2 then
 if maxim-minim=1 then LRT(2,`There is at most one more generator besides the`, minim,`known; conjecturally,`);LRT(2,`there is exactly one more.`);fi:
 if maxim-minim>1 then LRT(2,`The rank is between`,minim,`and`,maxim);fi:
else
 if maxim=1 then LRT(2,`There is at most one generator; conjecturally, there is exactly one.`);fi:
 if maxim>1 then LRT(2,`The rank is at most`,maxim,`and conjecturally of that parity.`);fi:
fi:
RETURN(maxim):
#compare to  nops(W)  (nops(W)/2 if TORSION=2):  what did we gain?
end:

ONED1:=proc(d1, A, d1d2) #extra args = seed to select other points on conic
local s,AA,dd:
#
#Check to see if a locally solvable homog space for an EC w/ 2-torsion
#is really solvable, at least as far as 2nd 2-descent
#
if nargs=0 then
	lprint(`ONED1(d1,A,B): Check to see whether the one homogeneous space`);
	lprint(`y^2=d1x^4+Ax^2+(B/d1) of E:y^2=x(x^2+Ax+B) admits higher descent.`);
	lprint(`An additional rational argument can be supplied to change the reference point`);
	lprint(`on the associated conic, thus yielding different factoring challenges, e.g.`);
	lprint(`For more info try  HELP(): `); RETURN():
fi:
LRT(4,`You have asked for a decision on the`,d1,` -homogeneous space for parameters`):
LRT(4, d1d2,` and `,A):
PRETTIFY(A,d1d2):AA:="[1]:dd:=""[2]:
if TORSION=0 then RETURN({}):fi:
if not(type(d1,integer) and type(dd/d1,integer)) then
	ERR(d1,`is not an integer divisor of`,dd); RETURN({}):fi:
if nargs>3 then RETURN(GUTS(d1,AA,dd/d1,args[4])): 
	else RETURN(GUTS(d1,AA,dd/d1)):fi
#maybe loop around looking for  args[4]  which give solvable factoriz. probls.?
end:

MANYD1:=proc(S,A,d1d2) 
local T,s,n,W,AA,dd:
#
# Try higher descent on many homog. spaces associated to a single  EC.  
#
if nargs=0 then
	lprint(`MANYD1(S,A,B): Check to see which of the homogeneous spaces`);
	lprint(`y^2=d1x^4+Ax^2+(B/d1) of E:y^2=x(x^2+Ax+B) admit higher descent,`);
	lprint(`as  d1  ranges over  S.`);
	lprint(`For more info try  HELP(): `); RETURN():
fi:
PRETTIFY(A,d1d2):AA:="[1]:dd:=""[2]:if TORSION=0 then RETURN({}):fi:
W:={}:
T:=sort([op(S)]):#I like to keep order uniform
for s in T do 
    if not(type(s,integer) and type(d1d2/s,integer)) then
	ERR(s,`is not an integer divisor of`,d1d2); RETURN({}):fi:
    n:=GUTS(s,A,d1d2/s):if n>0 then W:={op(W),s}:fi:od:
LRT(4,` _,-'"``-._,-'"``-._,-'"``-._,-'"``-._,-'"``-._,-'"``-._,-'"``-._,-'"``-._`);
if nops(S)=nops(W) then 
 LRT(4,`All`,nops(S), `original homogeneous spaces may contain rational points.`); 
else
 LRT(4,`Of the`,nops(S),`original homogeneous spaces presented,`,nops(W),`of them`);
 LRT(4,`may contain rational points, namely`,W);fi:
RETURN(W);
end:

#########################################################################
# In the routines which follow, I will assume various global parameters #
# have been appropriately set, and that data have been pre-screened for #
# various conditions (usually, for example, inputs are now required to  #
# be integral -- not required in front end). Use with caution.          #
#########################################################################

#CONTENTS:
#  OUTPUT FORMATTING
#  USUAL (FIRST) TWO-DESCENT 
#  GUTS OF 2nd 2-DESCENT
#  SWIFT CONIC CERTIFICATION uses:
#     LINEAR ALGEBRA MOD 2
#  POINTS ON CONICS  uses:
#     NUMBER THEORY BASICS (e.g. msqrt)
#  QUARTICS - LOCAL SOLVABILITY, REDUCTION, POINT SEARCH
#  NEXT-STEP DESCENT (PAIRS OF CONICS) - same themes? uses:
#     GAUSSIAN ARITHMETIC
#  INTEGER FACTORING ROUTINES
#  HELP FILES 

############################ OUTPUT FORMATTING ############################

PRT:=proc() local i:
if yaklevel>=args[1] then print(seq(args[i],i=2..nargs)); fi: end:
 #-1=totally silent (except HELP) 
 #0=error messages only 
 #2=Front End's output, Errors,  and Warnings
 #3=add timing information
 #4=add basic description of operation (mostly in GUTS)
 #5=extra data about computations (e.g. prepare for searches on new homog.sp.)
 #6=really boring data (good for debugging)
 #8=only adds witnessing primes for non-locally-solvable conics
 #10=shows entrance to every major procedure

LRT:=proc() local i:
if yaklevel>=args[1] then lprint(seq(args[i],i=2..nargs)); fi: end:

ERR:=proc() local i:
if yaklevel>=0 then print(`Error: `,seq(args[i],i=1..nargs)):fi:end:

WARN:=proc() local i:
if yaklevel>=2 then lprint(`Warning: `,seq(args[i],i=1..nargs)):fi:end:

PREPEXTERNAL:=proc() #prepare a file for an external factoring program.
local i:
writeto(`data`);
lprint(`1`);#I happen to use the LIDIA factorizer now.
for i in FATALSET do lprint(i);od;
lprint(`0`);
writeto(terminal):
end:

ADDREQUEST:=proc()
global SEARCHES:
if SEARCHES=0 then 
	writeto(`ratpoint.in`);
	lprint(1): #"verbose?"
else appendto(`ratpoint.in`):fi:
lprint(args);#"coefficients?"
lprint(0);#"roots known?"
lprint(8);#"search limit?"
writeto(terminal):
SEARCHES:=SEARCHES+1:
end:

END:=proc()
local K,i:
global KNOWNPR:
K:=KNOWNPR:
read(`prims.m`):#I might have other concurrent sessions updating this file too!
KNOWNPR:=KNOWNPR union K:
save(KNOWNPR, `prims.m`):
LRT(4,`Saving`,nops(KNOWNPR),`primes`);
LRT(3,`You spent`,HASSETIME,`seconds finding divisors allowing solvable conics.`);
LRT(3,`You spent`,NORMTIME,`seconds finding points on those conics.`):
LRT(3,`You spent`,LSTIME,`seconds deciding local solvability a la Siksek.`):
LRT(3,`You spent`,QTIME,`seconds improving the locally solvable quartics.`):
LRT(3,`You spent`,FACTTIME,`seconds factoring inside IFACTORS.`):
LRT(3,`(The last category overlaps the others.)`);
if factorproc=`external` and nops(FATALSET)>0 then
 LRT(4,`Here are the composite numbers you really need to factor:`);
 for i in FATALSET do LRT(4,i);od;
 PREPEXTERNAL():
 fi;
if WANTSEARCH and SEARCHES>0 then
 appendto(`ratpoint.in`);
 lprint(0,0,0,0,0):
 writeto(terminal):
fi:
quit:
end:

####################  USUAL (FIRST) TWO-DESCENT ##############################

PRETTIFY:=proc(A,d1d2)
local dd,AA,f,i:
global TORSION:
# Establish some parameters, fix data types to integers.
# (Could do more here, e.g. allow y^2=cubic with root <> 0)
LRT(10,`Entering prettify`,args):
if dd=0 then ERR(`Curve is singular`):TORSION:=0:RETURN([0,0]):fi:
dd:=numer(d1d2):
f:=IFACTORS(denom(d1d2))[2]:for i in f do dd:=dd*i[1]^(3-modp(i[2]-1,4)):od:
#so now it's an integer, and differs from orig. by a fourth power.
f:=sqrt(dd/d1d2):AA:=f*A:
if not type(AA,integer) then dd:=denom(AA)^4*dd:AA:=AA*denom(AA)^2:fi:
f:=AA^2-4*dd:
if f=0 then ERR(`Curve is singular`):TORSION:=0:RETURN([0,0]):fi:
if issqr(f) then TORSION:=2:
	LRT(4,`Curve contains 3 elements of order 2`);
	WARN(`Curves with such 2-torsion known to cause front-end problems.`):
else TORSION:=1:fi:
#Should probably work harder to make curve minimal. Consider  SQRPART(AA)...
RETURN([AA,dd]):
end:

FIND_D1S:=proc(A,d1d2)
local S,T,W,ct,d3,p,d1,S2,S1,S0,rk, D1SET,ht:
global LIKELIES,HASSETIME:
#
# Look for divisors  d1 of d1d2   for which quartic y^2=d1x^4+Ax^2+d2 is
# everywhere locally solvable.
#	NOTE: something is wrong when TORSION > 1 -- picks up useless d1's.
#
LRT(10,`Entering find_d1s`,args);
LRT(4,`First checking for divisors  d1  which make the quadratic solvable.`);
ht:=time():
S0:=FACINT(d1d2): CHKSNAG():
S0:=S0 minus {max(op(S0))}:#we can avoid one prime, knowing there's torsion
S1:=ODDINT(A^2-4*d1d2): CHKSNAG(): #S1={} means Z_2^2  < EC.
S2:=FACINT(A^2-4*d1d2):
if nops(S0) < 6 then #just an arbitrary cutoff to improve speed
 T:=sort([op(combinat[powerset](S0 union {-1}))]):
 W:={}:ct:=0:
 for S in T do 
  p:=WITNESS(S1,S):
  if p=1 then PRT(6,S,` : conic impossible to solve in reals`); 
  elif p>0 then PRT(6,S,` : conic impossible at p=`,p); 
  else W:={op(W),S}:
  fi:
  ct:=ct+1:if (ct mod 256)=0 then LRT(4,`(`,ct,`divisors checked)`);fi;
 od:
else
 S:=WHOSANORM(S1,S0 union {-1}):W:=FOLDUP(vector(nops(S[2]),0),S[1],S[2],2):
fi:
HASSETIME:=HASSETIME+time()-ht:
LRT(5,nops(W),`divisors d1 allow the quadratic equation to be solved. Of these,`);
D1SET:={}:LIKELIES:=[2]:
T:=S0 union S2 union {2}:#"Everywhere" means these p.
for S in sort([op(W)]) do
 d1:=1:for p in S do d1:=d1*p:od:
 if LOCALLYSOLV(d1,0,A,0,d1d2/d1,T) then D1SET:=D1SET union {d1}:
    PRT(5,`d1=`,d1,`gives a quartic solvable everywhere`):
    if SNAG >0 then PRT(5,`except perhaps at primes in`,SNAG);
       PRT(5,`Higher descent computations will fail without this data.`):fi:
    fi:
od:
D1SET:=D1SET minus {1};
LRT(4,`The following`,nops(D1SET),`divisors  d1  allow points on the homogeneous space:`);
PRT(4,D1SET);
rk:=round(evalf(log(nops(D1SET)+1)/log(2))):
LRT(4,`Thus the rank of the (E/tors)/im(phi) is at most`,rk);
LRT(4,`(and conjecturally of that parity).`);
RETURN(D1SET);
end:

######################## GUTS OF 2nd 2-DESCENT ##############################

GUTS:=proc(d1, A, d2)
local bigd,sol,z0,x0,y0,S0,S1,Sa,S2,S3,Sx,Sz,Sxz,d3s,qr,num,ht:
global FATAL,BADSET,HASSETIME:
#
# Main procedure: check to see if an "even"  homog. space can be re-descended.
# Note: different choices of an extra argument yield different factoring probs.
#
LRT(10,`Entering guts`,args);
LRT(4,` =================================================================== `);
LRT(4,`Attempting second 2-descent on the homogeneous space`);
PRT(4,`y^2 = `,d1,`x^4 + `,A,`x^2 + `,d2);
if issqr(d2) then 
  LRT(4,`Clearly`,[0,sqrt(d2)],`is a point on this space.`);
  RETURN(1):fi:
#Failure to trap these can allow  NEWHOMOG to create a _cubic_, not quartic.
bigd:=A^2-4*d1*d2:
FATAL:=false:
S1:=ODDINT(d1):CHKSNAG():
S2:=ODDINT(d2):CHKSNAG():
S3:=FACINT(bigd):CHKSNAG():
if FATAL then RETURN(0):fi: #That's _really_ Fatal! So assume FATAL is false..
BADSET:=FACINT(d1) union FACINT(d2) union S3:
LRT(4,`Given such a point, X=x^2 and Y=y will satisfy`);
PRT(4,Y^2 = d1*X^2 + A*X + d2);
LRT(4,`One point on this conic has [X,Y]=`);
#
sol:=NORM(d1,bigd,A);
if FATAL then ERR(`Unable to find a point on conic!`);RETURN(0):fi:
[(sol[1]-A*sol[3])/2/d1,sol[2]/2,sol[3]]:
sol:=["[1]*denom("[1]),"[2]*denom("[1]),"[3]*denom("[1])]:
x0:=sol[1]:y0:=sol[2]:z0:=sol[3]:
sol:=ilcm(denom(x0),denom(y0),denom(z0)):x0:=x0*sol:y0:=y0*sol:z0:=z0*sol:
#
if z0=0 then 
  ERR(`Homogeneous space seems not everywhere locally solvable.`):RETURN(0):fi:
if issqr(x0*z0) then
 LRT(4,`What luck;`,[sqrt(x0/z0),y0/z0],`is a point on the homogeneous space`);
 RETURN(1):#a positive RETURN means there's a possible point
fi: #in particular we may assume now that  x0<>0, z0<>0.
if nargs>3 then qr:=args[4]: #try a new spot on the curve 
  x0:=((x0/z0)*qr^2-2*(y0/z0)*qr+(A+d1*x0/z0))/(qr^2-d1):
  y0:=sqrt(subs(x=x0,d1*x^2+A*x+d2)):
  z0:=lcm(denom(x0),denom(y0)):
  x0:=x0*z0:y0:=y0*z0:
  fi:
PRT(4,[x0/z0,y0/z0]);
LRT(4,`All other points then have  X  given by  Xn/Xd  where (for some  t)`);
LRT(4,` `):
PRT(4,Xn=(x0/z0)*t^2-2*(y0/z0)*t+(A+d1*x0/z0));
PRT(4,Xd=t^2-d1);
LRT(4,`The homogeneous space has a point iff we may select  t  to make  Xn/Xd  a`); 
LRT(4,`rational square. This requires we select  t  and a square-free integer d3`);
LRT(4,`so that both Xn=d3*square and Xd=d3*square. Clearing denominators shows`);
LRT(4,`d3  will be an integer factor of the resultant of  Xn  and  Xd, hence`);
LRT(4,`a product of primes from`);
Sx:=ODDINT(x0):CHKSNAG():
Sz:=ODDINT(z0):CHKSNAG():
if FATAL then ERR(`Pass GUTS a 4th argument to circumvent factoring trouble?`);
	RETURN(0);fi:
BADSET:=BADSET union FACINT(z0):
Sxz:=symmdiff(Sx,Sz):
S0:=S3 union {-1} union FACINT(z0):#OK primes for  d3 -- divisors of resultant
LRT(4,S0):
#
ht:=time():
LRT(4,`The d3 for which such a  t  exists are limited to these:`);
if nops(S0)< 6 then #    the cutoff "6" is arbitrary; adjust as you think best.
 Sx:= TWOCONICS1(S1, S2, Sxz, S0);#success in both conics separately
 #But TWOCONICS1 fails if S0>16 (Maple can't index a set with >2^16 elements)
else
 Sx:= TWOCONICS2(S1, S2, Sxz, S0);# faster when  S0  is large
fi: 
#if nops(S0) is large, maybe add a parameter to ONED1 & see if S0 get smaller?
if nops(Sx) >0 then
 LRT(4,nops(Sx),`divisors d3 making both conics separately solvable.`);
else LRT(4,`none also allow the numerator to be solved.`);fi;
HASSETIME:=HASSETIME+time()-ht:
Sx:=SIMULT_CONIC(Sx, [d1,A,x0,y0,z0]); #Which d3 are _simultaneously_ loc.solv.
num:=AMBIG+nops(Sx):
#
if num>0 then 
# LRT(4,`We have found:`);
# if AMBIG>0 then  
#  LRT(4,`    `,AMBIG,` divisors which could not be checked (factoring trouble)`);fi:
# if nops(Sx)>0 then 
#  LRT(4,`    `,nops(Sx)-AMBIG,` definite descendant curves.`);fi:
#--unnecessary; the true value of nops(Sx) is zero or #{locallyOK d1s for E'}.
 LRT(4,`The original homogeneous space may have rational points.`);
else LRT(4,`The original homogeneous space has no rational points.`);fi;
RETURN(num):
end:

TWOCONICS1:=proc(S1,S2,Sxz,S3)
local S,T,W,ct,d3,p:
#
# Look for divisors  d3  in  S3  for which both (Xd=) t^2-d1 =d3*square and 
# (Xn=) (x0/z0) [(t-y0/x0)^2 - d2 (z0/x0)^2]=d3*square are solvable
#
LRT(10,`Entering twoconics1`);
LRT(5,`First checking divisors  d3  which make both quadratics separately solvable.`);
T:=combinat[powerset](S3):#should say  T={{}} if S3={}...
W:={}:ct:=0:
for S in T do 
 p:=WITNESS(S1,S):
 if p=1 then PRT(8,S,` : Xd impossible to solve in reals`); 
 elif p>0 then PRT(8,S,` : Xd impossible at p=`,p); 
 else W:=W union {S}:
 fi:
 ct:=ct+1:if (ct mod 256)=0 then LRT(4,`(`,ct,`divisors checked)`);fi;
od:
LRT(4,nops(W),`divisors d3 allow the denominator equation to be solved. Of these,`);
T:=W:W:={}:ct:=0:
for S in T do 
 p:=WITNESS(S2,symmdiff(Sxz,S)):
 if p=1 then PRT(8,S,` : Xn impossible to solve in reals`); 
 elif p>0 then PRT(8,S,` : Xn impossible at p=`,p); 
 else LRT(5,S,` : at least makes them separately solvable. `);
  d3:=1:for p in S do d3:=d3*p:od:
  W:=W union {d3}:
 fi:
 ct:=ct+1:if (ct mod 256)=0 then LRT(4,`(`,ct,`divisors checked)`);fi;
od:
RETURN(W):
end:

#TWOCONICS2, a substitute for large  S3, is discussed below

SIMULT_CONIC:=proc(S,pt)
local W,T,ct,d3,p,eq5,ee,A0,t0,z0:
global SNAG,FATAL,AMBIG,LIKELIES;
#
# Look for divisors  d3  in  S3  for which (Xd=) t^2-d1=d3*square and
# (Xn=) (x0/z0) [(t-y0/x0)^2 - d2 (z0/x0)^2] =d3*square  have a common
# solution  t  at all primes.
#
LRT(10,`Entering simult_conic`);
T:=sort([op(S)]):#guarantee uniformity of runs
LIKELIES:=[2]:AMBIG:=0:W:={}:ct:=0:
for d3 in T do
  ee:=NEWHOMOG(d3,op(pt)); #adding an extra parameter might help factor disc...
  if FATAL then
   LRT(4,`Warning: d3=`,d3);
   LRT(4,`leads to solvable conics but I cannot form quartic without those factors.`);
   FATAL:=false:
   AMBIG:=AMBIG+1:
  else
   LRT(5,`Solutions to making  Xd=d3*square are parameterized by`);
   LRT(5,` `);
   PRT(5,t=(d3*ee[1]*u^2-2*d3*ee[2]*u+ee[1])/(d3*u^2-1));
   LRT(5,`We substitute into  Xn  and see that now this must be square:`);
   eq5:=ee[3]:
   PRT(5,eq5[1]*u^4+eq5[2]*u^3+eq5[3]*u^2+eq5[4]*u+eq5[5]);
   if LOCALLYSOLV(op(eq5),BADSET union FACINT(denom(ee[1]))) then
    W:=W union {d3}:
    PRT(4,`d3=`,d3,`gives a quartic solvable everywhere`):
    if SNAG>0 then 
     #This should never be necessary now that I know how discrim(eq5) factors.
     LRT(4,`(except perhaps at divisors of`,SNAG,`)`):
     AMBIG:=AMBIG+1:
    else
      t0:=ee[1]:A0:=ee[2]:z0:=pt[5]:
      LRT(4,`The quartic is`);
#      PRT(4,v^2=eq5[1]*u^4+eq5[2]*u^3+eq5[3]*u^2+eq5[4]*u+eq5[5]);
#      LRT(4,`If [u,v] is a point on this curve, then `);
#      PRT(4, X=simplify(v/(d3*z0*denom(t0)*(A0*d3*u^2-2*t0*u+A0))));
#      LRT(4,`gives a point on the original homogeneous space.`);
#      LRT(4,`A simpler, equivalent quartic is`);
      LRT(4,IJREDUCE(op(eq5)) );
      if JUSTONE then break fi: 
          #i.e. if this d3 is really loc. solv., skip others.
          #ANy d3 works in principle _if_ Sha=0 on 2-isog. curve; but in
          #any case, some d3's will give quartics with smaller points
    fi:
   fi:
  fi:
  ct:=ct+1: if modp(ct,16)=0 then LRT(4,`(`,ct,`divisors checked)`):fi:
od:
RETURN(W);
end:

NEWHOMOG:=proc(d3,d1,A,x0,y0,z0)
local u,v,w,t,eq5,eq5c,ff,i,p,t0,A0,got,m1;
global FATAL:
#
#Calculate  new quartic for a d3 with both conics solvable.
#Use extra  parameter  m1 to get different points (possibly affects factoriz)
#
LRT(10,`Entering newhomog`);
FATAL:=false:
got:=NORM(d1,d3):if FATAL then RETURN([0,0,[0,0,0,0,0]]):fi:
if got[2]=0 then #can only happen when  d3=1!
 got:=[(d1+1),2,(d1-1)]:fi:
t0:=got[1]/got[2];A0:=got[3]/got[2]:#so t0^2-d1=d3*A0^2
if nargs>6 then m1:=args[7]:
  t0:=(t0*d3*m1^2-2*d3*A0*m1+t0)/(d3*m1^2-1):
  A0:=sqrt((t0^2-d1)/d3):
  fi:
t:=(t0*d3*nu^2-2*d3*A0*nu+t0)/(d3*nu^2-1): #makes Xd= t^2-1=d3*V^2, where V=
#(nu^2*A0*d3-2*nu*t0+A0)/(d3*nu^2-1)
eq5:=simplify((x0*t^2-2*y0*t+d1*x0+A*z0)*d3*z0);#choose  nu  to make square?
#This is d3*z0^2*Xn, so if=U^2 then Xn/Xd=(U^2/d3/z0^2)/(d3*V^2)=(U/d3/z0/V)^2
eq5:=expand(simplify(eq5*(d3*nu^2-1)^2)):#"expand" cuz could be ()*().[BSD]
#If this is U'^2, U=U'/(d3nu^2-1), Xn/Xd=(U'/d3/z0/numer(V))^2
# eq5 should now be a _quartic polynomial_  in  nu. Note: x0*z0*d3<>0.
#The only denominators are those of t0 and A0 (common), squared:
eq5:=eq5*denom(t0)^2;#now it should be _integral_, too.
#If this is U"^2, then (Xn/Xd)=[U"/(d3*z0*denom(t0)*numer(V))]^2.
eq5c:=[coeff(eq5,nu,4),coeff(eq5,nu,3),coeff(eq5,nu,2),coeff(eq5,nu,1),coeff(eq5,nu,0)];
#This is what must be square (upon suitable choice of  nu) 
#The following should never occur (but would hang  LOCALLYSOLV):
if eq5c[1]=0 then 
  ERR(`Newhomog produced a nonquartic`):RETURN([0,0,[0,0,0,0,0]]):fi:
RETURN([t0,A0,eq5c]);
end:

################### SWIFT CONIC CERTIFICATION ###########################
#Faster way to determine divisors which allow quadratics to be solved.
#Still could use some better linear-algebra analysis; see below.

WITNESS:=proc(S1,S2)
local p,k,d,h, S1a, S2a, S3:
#
#Decide if u^2-dv^2=hw^2  has rational solutions  u,v,w.
#input: sets S1, S2 of primes (or -1) whose products are d, h resp. mod squares
#output: p=0 => solvable; p=1 => can't solve in  R;  p>1: can't solve in Q_p.
#
LRT(10,`Entering witness`);
S3:=S1 intersect S2:
S1a:=S1 minus S3:
S2a:=S2 minus S3:
d:=1:for p in S1a do d:=d*p: od:
h:=1:for p in S2a do h:=h*p: od:
k:=-1:for p in S3 do k:=k*p: od:
#We're really solving  kx^2+dy^2+hz^2=0.
if (k<0) and (h<0) and (d<0) then RETURN(1):fi: #can't really happen
if (k>0) and (h>0) and (d>0) then RETURN(1):fi: #"-1" in both S1,S2
for p in S1a do if L(-k,p)=-L(h,p) then RETURN(p):fi:od:
for p in S2a do if L(-d,p)=-L(k,p) then RETURN(p):fi:od:
for p in S3  do if L(-h,p)=-L(d,p) then RETURN(p):fi:od:
#Note: this test is _incorrect_ mod 2, but an equ. can't be locally solvable
#_everywhere_ else and yet not at 2.
RETURN(0):
end:

WHOSANORM:=proc(S1,S3)
local d1,d3,p,S3a,S3b,S0,S1b,T0,T1,T3,t0,r,c,R,C,MN,MD,Basis,i,j,v,u,U,Q1,Q2,B,v1,v2,w,LFT,RGHT:
#
# Look for divisors  d3  in  S3  for which  t^2-d1=d3*square  is solvable,
# where  d1  is the product of the primes (or (-1)) in S1.
#
# I must be missing something simple here: the set of such d3 is a vector
# subspace of  Q/Q^2, (it's the image of Q[sqrt(d1)] under the norm map)
# but I can't find a way to describe it all as a kernel
# or something else easily computed in linear algebra. The following is
# close (it's inefficient only to the extent that  S1 intersect S3 > {})
#
LRT(10,`Entering whosanorm`);
if member(-1,S1) then S3a:=S3 minus {-1}: else S3a:=S3:fi: 
#Ensures real solvability; now make solvable mod all odd  p  dividing d1*d3.
S0:=S1 intersect S3a:S1b:=S1 minus S0:
#Can only use primes in  S3a  mod which  d1 is a square:
d1:=1:for p in S1 do d1:=d1*p:od:
S3b:={}:for p in S3a minus S0 do if L(d1,p)=+1 then S3b:=S3b union {p}: fi:od:
#Even if -1 is in S1a, it won't affect further choices of d3's, so do
S1b:=S1b minus {-1}:
#Three rel. prime sets of primes. Fix an ordering of each:
T0:=[op(S0)]:T0:=sort(T0):t0:=nops(T0);
T1:=[op(S1b)]:T1:=sort(T1):R:=[op(T0),op(T1)]:r:=nops(R):#rows<->primes of d1
T3:=[op(S3b)]:T3:=sort(T3):C:=[op(T0),op(T3)]:c:=nops(C):#cols<->primes of d3
if r=0 then #r=0 happens if  d1=-1.
 Basis:={}:
 for i from 1 to c do v:=vector(c,0):v[i]:=1:
  Basis:=Basis union {vector([seq(v[j],j=1..c)])}:od:
 RETURN([Basis,C]):
else #i.e., if r>0 
 if c>0 then MN:=matrix(r,c):
  for i from 1 to r do for j from 1 to c do
    MN[i,j]:=(1-L(C[j],R[i]))/2:#test the prime divs. of N=d3: =sq. mod (p|d1)?
  od:od:fi: #Jacobi symbols as a 2x2-partitioned matrix over field  F_2
 if t0>0 then LFT:=matrix(r,t0):
  for i from 1 to r do for j from 1 to t0 do
    LFT[i,j]:=MN[i,j]:
  od:od: fi:
 if c>t0 then RGHT:=matrix(r,c-t0):
  for i from 1 to r do for j from 1 to c-t0 do
    RGHT[i,j]:=MN[i,j+t0]:
  od:od: fi:
fi:
if t0>0 then
 MD:=vector(t0):for j from 1 to t0 do MD[j]:=(1-L(-d1/C[j],C[j]))/2:od:fi:
#(Additively we could write L(x/p,p)"="L(x,p)-L(p,p)=L(x,p), but L isnt really
#additive: L(x,p)=0 for p|x since L=0 just means x is a square.)
#
#Now want the (Z/2-vector space!) of column vectors  v = (v1,v2) such that:
#"L(D,N)=+1": we already ensured this by limiting S3b
#"L(N,D)=+1": matrix product  MN*v=(*,0) where also "*"[i]=0 whenever v1[i]=0
#"L(-D/gcd * N/gcd , gcd)=+1": expressed with  LFT and RGHT (see below)
#We record answer via a basis: first all those in (0,*), then a preimage of 
#image of one in (*,0).
#
#First find all solutions with g=gcd(d1,d3)=1, i.e, v = (0,*).
#In that case, just need  L(N,p)=+1 for all  p | d1: 
if (c-t0=0 or r=0) then U:={} else U:=NULLSPACE(RGHT): fi:
Basis:={}:for v in U do 
  Basis:=Basis union {vector([seq(0,i=1..t0),seq(v[i],i=1..c-t0)])};od:
#
#Now complete to a basis by looking for vectors v=(v1,v2) for each v1<>0.
#For p | (d1/g) need L(N/g,p)=L(g,p) so L(N,p)=0
#    (In particular, v  will lie in the nullspace of the (r-t0) x c "bottom"
#    (of  MN. So we could use more lin. algebra to speed this up. Not done yet.
#For p | g need L(N/g,p)=L(-d1/g,p), so L(N/p,p)=L(-d1/p,p)
#In either case, the LHS is the sum of  L(q,p) for all  q|N
#That is, we want  MN*v to be the vector whose ith coord is  MD[i] unless
#the ith prime divides  d1/g, ie isnt in g ie isnt in d3 ie v1_i=0.
#Now consider all possible  gcd's  g<>1:
if t0>0 then U:=powerset(t0) minus {{}}; else U:={}:fi:
for u in  U  do  #u  will pick out the set of primes in  g = gcd(d1,d3).
  B:=vector(r,0): for i in u do B[i]:=MD[i]:od: #so that's the target
  #want  MN*v=target but want v=(u,v2), so want LFT*v1+RGHT*v2=target:
  v1:=vector(t0,0): for i in u do v1[i]:=1:od:
  Q1:=multiply(LFT,v1);Q2:=[seq(B[i]-Q1[i] mod 2,i=1..r)]:
  if (r=0 or t0=c) then v2:=LSOLV2(Q2) else
  v2:=LSOLV(RGHT,Q2)[1]:fi: #just need _one_ solution giving correct gcd(d1,N).
  if nops(v2)>0 then  #i.e. if there IS a solution,
    v:=vector(c,0):
    for j from 1 to t0 do v[j]:=v1[j]:od:
    for j from 1 to c-t0 do v[j+t0]:=v2[j]: od:
    Basis:=Basis union {vector([seq(v[i],i=1..c)])}:
  fi:
od:
if nops(Basis)>0 then Basis:=TRIM(Basis): fi:
RETURN([Basis,C]);
end:

FOLDUP:=proc(partic,bas,primes,flag)
local d,i,U,u,V,v:
# Convert vectors "d3" to sets of primes or products of primes.
#
LRT(10,`Entering foldup`);
V:={}:if nops(bas)=0 then U:={{}} else U:=powerset(nops(bas)):fi:
for u in U do
 v:=partic:for i in u do v:=ADDVECT(v,bas[i]):od:
 V:=V union {vector([seq(v[i],i=1..nops(primes))])}:
od:
U:={}:for v in V do
 if flag=1 then
  d:=1:for i from 1 to nops(primes) do if v[i]=1 then d:=d*primes[i]:fi:od:fi:
 if flag=2 then
  d:={}:for i  from 1 to nops(primes) 
      do if v[i]=1 then d:={op(d),primes[i]}:fi:od:fi:
 U:={op(U),d}:
od:
RETURN(U):
end:

#Now use this procedure to make a fast plug-in replacement for  TWOCONICS1:

TWOCONICS2:=proc(S1,S2,Sxz,S3)
local S,T,W,d3,Z,Ordr1,Bas1,Ordr2,Bas2,n1,n2,n,Ordr,cnvrt1,cnvrt2,i,j,Bas1a,Bas2a,Bas4,Bas5,u,v,w,target,Partic,M,U,u,n4,n6,Bas6:
#
# Look for divisors  d3  in  S3  for which both (Xd=) t^2-d1 =d3*square and 
# (Xn=) (x0/z0) [(t-y0/x0)^2 - d2 (z0/x0)^2]=d3*square are solvable
#
LRT(10,`Entering twoconics2`);
LRT(5,`First checking divisors  d3  which make both quadratics separately solvable.`);
Z:=WHOSANORM(S1,S3):
Bas1:=Z[1]:Ordr1:=Z[2]:
n1:=nops(Bas1):
LRT(4,2^n1,`divisors d3 allow the denominator equation to be solved. Of these,`);
Z:=WHOSANORM(S2,Sxz union S3);
Bas2:=Z[1]:Ordr2:=Z[2]:
n2:=nops(Bas2):
#
#Bring the two bases into alignment first:
Ordr:=sort([op({op(Ordr1)} union {op(Ordr2)})]):n:=nops(Ordr):
cnvrt1:=vector(nops(Ordr1)):
for i from 1 to nops(Ordr1) do
 for j from 1 to n do if Ordr[j]=Ordr1[i] then cnvrt1[i]:=j:fi:od:od:
cnvrt2:=vector(nops(Ordr2)):
for i from 1 to nops(Ordr2) do
 for j from 1 to n do if Ordr[j]=Ordr2[i] then cnvrt2[i]:=j:fi:od:od:
Bas1a:={}:for v in Bas1 do
 w:=vector(n,0):for i from 1 to nops(Ordr1) do w[cnvrt1[i]]:=v[i]:od:
 Bas1a:=Bas1a union {vector([seq(w[i],i=1..n)])}:od:
Bas2a:={}:for v in Bas2 do
 w:=vector(n,0):for i from 1 to nops(Ordr2) do w[cnvrt2[i]]:=v[i]:od:
 Bas2a:=Bas2a union {vector([seq(w[i],i=1..n)])}:od:
#
#Now we want to know which vectors in (span Bas2a)n(span S3) have a sum with 
# xz which is still in that first span, that is, need projection to
#(Sxz union S3)/(S3) to be exactly  xz.
target:=vector(n,0):
for i from 1 to n do if member(Ordr[i],Sxz) then target[i]:=1 fi:od:
M:=matrix(n, n2):
for j from 1 to n2 do for i from 1 to n do M[i,j]:=Bas2a[j][i]:od:od:
for i from 1 to n do 
 if member(Ordr[i],S3) then 
  for j from 1 to n2 do M[i,j]:=0:od:#OK, so this is a lame solution
  target[i]:=0:
 fi:
od:
Z:=LSOLV(M,target):
w:=Z[1]:
if w=[] then 
 LRT(5,`The second conic isn't even separately solvable!`):RETURN({}):fi:
Partic:=vector(n,0):
for i from 1 to n do if member(Ordr[i],Sxz) then Partic[i]:=1 fi:od:
for i from 1 to n2 do if w[i]=1 then Partic:=ADDVECT(Partic, Bas2a[i]):fi:od:
Bas4:={}:for v in Z[2] do
 w:=vector(n,0):
 for i from 1 to n2 do if v[i]=1 then w:=ADDVECT(w, Bas2a[i]):fi:od:
 Bas4:=Bas4 union {vector([seq(w[i],i=1..n)])}:
od:
n4:=nops(Bas4):
#so now the divisors of  S3  making x^2-d2 y^2= xz d3 w^2 solvable are those
#of the form (span of Bas4)+Partic.
LRT(6,`Second conic has points for precisely`,2^n4,`d3s.`):
#
#Now need to intersect a subspace with a coset. Solve Sum(xi vi)=Sum(yi wi)+P
M:=matrix(n,n1+n4);
for j from 1 to n1 do for i from 1 to n do M[i,j]:=Bas1a[j][i]:od:od:
for j from 1 to n4 do for i from 1 to n do M[i,j+n1]:=Bas4[j][i]:od:od:
Z:=LSOLV(M,Partic):
w:=Z[1]:
if w=[] then LRT(5,`No d3 makes both conics solvable`):RETURN({}):fi:
Partic:=vector(n,0):
for i from 1 to n1 do if w[i]=1 then Partic:=ADDVECT(Partic, Bas1a[i]):fi:od:
Bas5:={}:for v in Z[2] do
 w:=vector(n,0):
 for i from 1 to n1 do if v[i]=1 then w:=ADDVECT(w, Bas1a[i]):fi:od:
 Bas5:=Bas5 union {vector([seq(w[i],i=1..n)])}:
od:
Bas6:=TRIM(Bas5):
n6:=nops(Bas6):
LRT(5,`I find`,2^n6,`d3s make both conics solvable separately.`):
RETURN(FOLDUP(Partic, Bas6, Ordr,1)):
end:

################ LINEAR ALGEBRA MOD 2 ##############################
#Needed only for previous section of code.
#Much Maple code recreated due to what appears to be a missing line
#(fails to find nullspace of r x k matrix of rank r if k > r !)
#Note: there've got to be better primitives for exporting/adding/testing...
#vectors and matrices! (Why, for example, can't I say things like
#Basis:=Basis union {v}:  to adjoin a vector to a basis?)

TRIM:=proc(V)
local A,B,C,N,pivots,i,j,rk,Basis,l: 
# Take a set  V  of vectors mod 2 and return a minimal spanning set.
#
LRT(10,`Entering trim`);
l:=vectdim(V[1]):A:=matrix(nops(V),l):
for i from 1 to nops(V) do for j from 1 to l do A[i,j]:=V[i][j]:od:od:
B:=transpose(A);
N:=NULLSPACE(B);
if nops(N)=0 then RETURN(V):fi:
l:=vectdim(N[1]):C:=matrix(nops(N),l):
for i from 1 to nops(N) do for j from 1 to l do C[i,j]:=N[i][j]:od:od:
C:=GAUSSELIM(C);
rk:=C[2]:C:=C[1]:
pivots:={}:
for i from 1 to rk do
 j:=1:while C[i,j]=0 do j:=j+1:od:
 pivots:={op(pivots),j};
od;
Basis:={};
for i from 1 to nops(V) do 
	if not member(i,pivots) then
	   Basis:=Basis union {vector([seq(V[i][j],j=1..vectdim(V[i]))])}:fi:
od:
RETURN(Basis);
end:

NULLSPACE:=proc(A) #return a spanning set for nullspace of a matrix  A
local M,N,l,m,n,i,j,p,r,s,t,u,v,nullity;
#small modif of nullspace(), `Copyright 1993 by Waterloo Maple Software`;
LRT(10,`Entering nullspace`):
    M := GAUSSJORD(A):
    r:=M[2]: M:=M[1]:
    l := linalg[rowdim](A);
    m := linalg[coldim](A);
    n := m-r;
    if n = 0 then RETURN({}) fi;
    if 0 < r then s := array(1 .. r) fi;
    u := array(1 .. n);
    v := array(1 .. n);
    j := 1;
    for i to l while j <= m do
        while j <= m and M[i,j] = 0 do
            u[j-i+1] := j; v[j-i+1] := i; j := j+1;
        od;
        if i <= r then s[i] := j fi;
        j := j+1
    od;
while j<=m do u[j-l]:=j: v[j-l]:=i: j:=j+1: od:#this line was missing!
    N := array(1 .. n);
    for i to n do
        t := array(1 .. m);
        for j to m do  t[j] := 0 od;
        t[u[i]] := 1;
        for j to v[i]-1 do  t[s[j]] := (M[j,u[i]])  od;
        N[i] := op(t)
    od;
    {seq(eval(N[i]),i = 1 .. n)}
end:

GAUSSJORD:=proc(A)
local B,n,m,i,j,c,r,s,rank;
#small modif of gaussjord(), `Copyright 1993 by Waterloo Maple Software`;
LRT(10,`Entering gaussjord`);
    B := GAUSSELIM(A):
    rank:=B[2]: B:=B[1]:
    n := linalg['rowdim'](A);
    m := linalg['coldim'](A);
    r := 1;
    for c to m while r <= n do
        for i from r to n while B[i,c] = 0 do   od;
        if n < i then next fi;
        B[r,c] := 1;
        for i to r-1 do
            if B[i,c] = 0 then next fi; #so B[i,c]=1
            for j from c+1 to m do
                B[i,j] := (B[i,j]-B[r,j]) mod 2
            od;
            B[i,c] := 0
        od;
        r := r+1
    od;
    [op(B),rank]
end:

GAUSSELIM:=proc(A)
local B,d,n,m,i,j,c,r,t,rank;
#small modif of gausselim() `Copyright 1993 by Waterloo Maple Software`;
LRT(10,`Entering gausselim`);
    n := linalg['rowdim'](A);
    m := linalg['coldim'](A);
    B := array(1 .. n,1 .. m);
    for i to n do  for j to m do  B[i,j] :=modp(A[i,j],2)   od od;
    r := 1;
    for c to m while r <= n do
        for i from r to n while B[i,c] = 0 do   od;
        if n < i then d := 0; next fi;
        if i <> r then
            for j from c to m do
                t := B[i,j]; B[i,j] := B[r,j]; B[r,j] := t
            od
        fi;
        for i from r+1 to n do
            if B[i,c] = 0 then next fi; #so B[i,c]=1
            for j from c+1 to m do
                B[i,j] :=(B[i,j]-B[r,j]) mod 2
            od;
            B[i,c] := 0
        od;
        r := r+1
    od;
    rank := r-1;
    [op(B),rank]
end:

LSOLV:=proc(A,B)
local i,j,r,c,v,w,ww,wx,V,W,M:
#I don't see a primitve for linsolve mod 2. We can use NULLSPACE (or msolve)
#output: one solution (or "[]" if none), and basis of null(A).
LRT(10,`Entering lsolv`);
r := linalg['rowdim'](A);
c := linalg['coldim'](A);
M:=matrix(r,c+1):
for i from 1 to r do 
 for j from 1 to c do M[i,j]:=A[i,j]:od:
 M[i,c+1]:=B[i]:
od:
W:=NULLSPACE(M):
ww:=[]:V:={}:
for w in W do 
 if w[c+1]>0 then 
  if ww=[] then ww:=[seq(w[i],i=1..c)]: 
  else wx:=[seq(w[i]-ww[i] mod 2,i=1..c)] :V:=V union {op(wx)}:fi:
 else V:=V union {[seq(w[i],i=1..c)]}:fi:
od:
RETURN([ww,V]);
end:

LSOLV2:=proc(B) #"solve" linear system  "[empty | B ]", i.e. see if  B=0
local t, i:
LRT(10,`Entering lsolv2`);
t:=true:for i in B do if i<>0 then t:=false:fi:od:
if t then RETURN([0,0]) else RETURN([]);fi:
end:

ADDVECT:=proc(u,v)
local i:
LRT(10,`Entering addvect`);
RETURN(vector([seq(u[i]+v[i] mod 2, i=1..vectdim(u))])):end:

############################ POINTS ON CONICS ################################
#Considerable effort was expended while programming trying to find ways of
#finding "small" rational points on conics "efficiently". A call to  NORM
#will carry this out. It points to  NORMLLL  which is fast and produces 
#small points. In case  NORMLLL fails, (which hasn't yet happened) it 
#points to  NORMXa which isn't really well made and might have to call the
#inefficient  NORM0. I have some variant procedures under construction but
#won't invest time in them until weaknesses in these routines are known.
#Certainly this can be made polynomial-time and not use factoring.
#
#I also supply half-hearted routines which try to find smaller points in
#some cases on a more general conic. If it fails it reverts to regular
#NORMLLL. (See description of extra arguments in  NORM)
#
#  Sample: solve x^2-123452853y^2=1502226171486750244 rationally.

NORM:=proc(d1,d0) 
local x,nt,g,M1,M2,M3,s1,s0,e1,e0,rt:
global NORMTIME,FATAL:
#
# Return  [x1,x2,x3]  with x1^2-d1x2^2=d0x3^2, or else set FATAL.
# (4th arg, if present, requests  x1 = A x3  mod  d1 - see GUTS.)
#
LRT(10,`Entering norm`,args);
nt:=time():
if nargs>2 then rt:=args[3]:else rt:=0:fi:#stipulate that x1=rt*x3 mod d1
LRT(5,`Finding a norm equal to`,d0,`in Z[ sqrt(`,d1,`) ]`);
s1:=SQRPART(d1):s0:=SQRPART(d0):e1:=d1/s1:e0:=d0/s0:s1:=sqrt(s1):s0:=sqrt(s0):
#should have  s1=1  in our applications.
g:=gcd(e1,e0):
FATAL:=false:
#NORMLLL was working well until I tried adding a 4th arg so I use caution:
if nargs=2 or g>1 then x:=NORMLLL(g,-e1/g,-e0/g) else 
 #sorry, this is definitely buggy when  g>1 (shouldnt be)
 if modp(rt,g)>0 then rt:=0:fi:#would need  x3=0 mod g//rt -- skip it!
 if gcd(s0,e1)>1 then rt:=0:fi:#not sure how to xfer extra condition. Too lazy.
 if rt=0 then 
  WARN(`Declining the opportunity for greater efficiency in  NORM.`);
  x:=NORMLLL(g,-e1/g,-e0/g):
 else
  rt:=rt/g/s0 mod (e1/g):
  x:=NORMLLL(g,-e1/g,-e0/g,rt):#if this fails, try without extra arg.:
  if x[1]=0 and x[2]=0 then x:=NORMLLL(g,-e1/g,-e0/g):fi:
 fi:
fi:
if x[1]=0 and x[2]=0 then #something's wrong with  NORMLLL (never yet seen) :
 WARN(`NORMLLL only got down to`,x[3],`; invoking  NORMXa instead`):
 x:=NORMXa(d1,d0):#backup procedure is better tested...
 NORMTIME:=NORMTIME+time()-nt:RETURN(x):
fi:
if (g*x[1]*s1*s0)^2-d1*(x[2]*s0)^2-d0*(x[3]*s1)^2<>0 then 
  ERR(`Bad NORM calculations.`):FATAL:=true:RETURN([0,0,0]):fi:
NORMTIME:=NORMTIME+time()-nt:
RETURN([g*x[1]*s1*s0,x[2]*s0,x[3]*s1]):
end:

NORMLLL:=proc(a,b,c) 
local x,t,rt:
#
# Find small(est) solution to  ax^2 + by^2 + cz^2 = 0 (with x/z=args[4] mod b)
#
LRT(10,`Entering normlll`,args):
if nargs>3 then rt:=args[4]:
	#stipulates  x=(rt)z mod b ; assumes  rt^2=-c/a mod b, so:
	if modp(a*rt^2+c,b)>0 then 
		WARN(`NORMLLL was passed a bad 3rd arg; ignoring it`):rt:=0:
	fi:
else rt:=0:fi:
if modp(a*b*c,2)=1 then
	if modp(a+b,4)=0 then x:=ODDABC(a,b,c,0,rt,0):
	elif modp(b+c,4)=0 then x:=ODDABC(b,c,a,rt,0,0):x:=[x[3],x[1],x[2]]:
	else x:=ODDABC(c,a,b,0,0,rt):x:=[x[2],x[3],x[1]]:
	fi:
else #abc is even
	if modp(c,2)=0 then x:=EVENABC(a,b,c,0,rt,0):
	elif modp(b,2)=0 then x:=EVENABC(c,a,b,0,0,rt):x:=[x[2],x[3],x[1]]:
	else x:=EVENABC(b,c,a,rt,0,0):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:
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:

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
#
LRT(10,`Entering scalelattice`,args):
Digits:=3*(length(a)+length(b)+length(c)+5):
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]:
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:

ODDABC:=proc(a,b,c) 
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
LRT(10,`Entering oddabc`,args);
g:=MSQRT(-a*b,c):if nargs>3 and args[6]<>0 then g:=args[6]:fi:g:=g/b mod c:
h:=MSQRT(-c*a,b):if nargs>3 and args[5]<>0 then h:=args[5]:fi:h:=h/a mod b:
i:=MSQRT(-b*c,a):if nargs>3 and args[4]<>0 then i:=args[4]:fi: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
WARN(`NORMLLL only hit`,g/4/a/b/c,`(not zero) with`,a,b,c,`; trying other v`);
#example: ODDABC(1,-3221,-1216435);
if a*v2[1]^2+b*v2[2]^2+c*v2[3]^2=0 then RETURN(v2):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*(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 
	WARN(`Still failed; here's a (non-small!) solution:`);
	#from Borevich+Shafarevich p. 64
	g:=v[1]*v[3]+2*b*v[2]:
	h:=v[2]*v[3]-2*a*v[1]:
	i:=v[3]^2+4*a*b:
	RETURN([g,h,i]):fi:
WARN(`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) 
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,c/2  are odd
g:=MSQRT(-a*b,c):if nargs>3 and args[6]<>0 then g:=args[6]:fi:g:=g/b mod c:
h:=MSQRT(-c*a,b):if nargs>3 and args[5]<>0 then h:=args[5]:fi:h:=h/a mod b:
i:=MSQRT(-b*c,a):if nargs>3 and args[4]<>0 then i:=args[4]:fi: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
WARN(`NORMLLL only hit`,g/4/a/b/c,`(not zero) with`,a,b,c,`; trying other v`);
#example: EVENABC(33,-77483,2):
if a*v2[1]^2+b*v2[2]^2+c*v2[3]^2=0 then RETURN(v2):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*(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 
	WARN(`Still failed; here's a (non-small!) solution:`);
	#from Borevich+Shafarevich p. 64
	g:=v[1]*v[3]+2*b*v[2]:
	h:=v[2]*v[3]-2*a*v[1]:
	i:=v[3]^2+4*a*b:
	RETURN([g,h,i]):fi:
WARN(`A small solution can be found using work of B.Vallee (not coded,sorry)`);
RETURN([0,0,g/a/b/c/4]);
end:

# This backup routine was stable but (since it calls  NORM0 at end) not fast:

NORMXa:=proc(d1,d0) 
local g, T1, T2, T3,S,S1,S2,S3,x0,y0,k,k1,d0a,d0b,d0c,d0d,d0e,aus;
#work-alike for NORM0; should be faster and give smaller numbers.
#d1 assumed to be a square-free integer.
d0a:=SQRPART(denom(d0)):d0b:=denom(d0)/d0a:d0c:=d0*d0b^2*d0a:# now integer
d0d:=SQRPART(d0c):d0e:=d0c/d0d:# now squarefree, too
aus:=NORMYa(d1,d0e):#       aus=[x,y,k] where norm(x,y)=k*d0, |k| <= |d0e|
x0:=aus[1]:y0:=aus[2]:k:=aus[3]:
g:=SQRPART(k):k1:=k/g:
if k1=1 then RETURN([x0,y0,sqrt(g*d0a/d0d)*d0b]);fi:
if abs(k1)<abs(d1) then S:=NORMXa(k1,d1): if FATAL then RETURN([0,0,0]):fi:
   S1:=S[1]:S2:=S[3]:S3:=S[2]/sqrt(g):
 elif abs(k1)<abs(d0) then S:=NORMXa(d1,k1): if FATAL then RETURN([0,0,0]):fi:
   S1:=S[1]:S2:=S[2]:S3:=S[3]/sqrt(g):
 else WARN(`NORM0 invoked`);S:=NORM0(d1,k1):if FATAL then RETURN([0,0,0]):fi:
    S1:=S[1]:S2:=S[2]:S3:=S[3]/sqrt(g):fi:
if FATAL then RETURN([0,0,0]);fi:
T1:=x0*S1-y0*S2*d1; #       Now divide the things with norm  k*d0 and k.
T2:=-x0*S2+y0*S1;
T3:=S3*k*sqrt(d0a)*d0b/sqrt(d0d):
g:=ilcm(denom(T1),denom(T2),denom(T3));# courtesy: make integral...
T1:=T1*g:T2:=T2*g:T3:=T3*g:
g:=igcd(T1,T2,T3): #               # ... and relatively prime.
T1:=T1/g:T2:=T2/g:T3:=T3/g:
if not (0 =T1^2-d1*T2^2-d0*T3^2) then ERR(`NORMXa BAD`,d1,d0):fi:
RETURN([T1,T2,T3]);
end:

NORMYa:=proc(d1,d0e) 
local g, T1, T2, T3, S, x0,y0,best,trial,x1,y1:
#d1,d0e assumed to be square-free integers. Return [x,y,k] with
# x^2-d1 y^2 = k d0e, k integral, k <= d0e in abs value.
LRT(10,`Entering normya`,args);
if d0e=1 then RETURN([1,0,1]):fi:
x0:=d0e:y0:=0:g:=gcd(d1,d0e):
S:=MSQRT(d1,d0e/g):if FATAL then 
	ERR(`I thought d0 was factored? d0=`,d0):RETURN([0,0,0]):fi:
x1:=g*S:y1:=g:
trial:=abs((x1^2-d1*y1^2)/d0e);
#so  d0e  divides (d0e + 0*zeta ) and (x1 + 1*zeta ) in Z[zeta], zeta=sqrt(d1).
#We know Euclidean norm is almost never good enough for gcd work, but
#it's usually helpful anyway. Minimize 'best':=norm(a+b zeta)/d0e  in the ideal
#( d0e,  x1+zeta) as follows
best:=trial+1;
while trial<best do
 best:=trial;
 g:=DIVALGa(x0,y0,x1,y1,d1);
 x0:=x1:y0:=y1:x1:=g[1]:y1:=g[2]:
 trial:=abs((x1^2-d1*y1^2)/d0e);
 if trial=0 then trial:=best+1:fi:
od;
best:=(x0^2-d1*y0^2)/d0e;#no abs now
RETURN([x0,y0,best]):
end:

DIVALGa:=proc(x0,y0,x1,y1,dd)
local a,b,c,d,e,f,g:
LRT(10,`Entering divalga`,args);
a:=x0*x1-dd*y0*y1;
b:=-x0*y1+y0*x1;
c:=x1^2-dd*y1^2;
d:=round(a/c):e:=round(b/c):
f:=x0-(x1*d+y1*e*dd);
g:=y0-(x1*e+y1*d);
RETURN([f,g]);
end:

# original idea, taken from Maple's  isolve.
# (factors too much, yields big answers; NORMXa uses in case of emergency)

NORM0:=proc(d1,d0)
local A,B,C,bc,S,d0a,d0b,d0c,d0d,d0e;
LRT(10,`Entering norm0`,args);
#
#Here's what's left of isolve if you strip the unused clauses out
#for aplications of the form  isolve(u^2-d1*v^2=d0*w^2,{n1,n2,n3});
#with either  d0=bigd or d0=d3. I,e, we need  NORM(u/w + v/w sqrt(d1))=d0.
#Note: d1  assumed to be squarefree integer! 
#Also, equation is assumed to be locally solvable!(see WITNESS)
#Hardly need this now, except  NORMXa  will only go so far.
#
LRT(5,`Seeking a norm of`,d0,`in Z[sqrt(`,d1,`)].`);
d0a:=SQRPART(denom(d0)):d0b:=denom(d0)/d0a:d0c:=d0*d0b^2*d0a:# now integer
d0d:=SQRPART(d0c):d0e:=d0c/d0d:# now squarefree
bc:=igcd(-d1,-d0e):
A:=bc:
B:=-d1/bc:
C:=-d0e/bc:
#if not ( L(-B*C,abs(A)) =1 and L(-A*B,abs(C))=1 and L(-C*A, abs(B))=1 and
#	not(d1<0 and d0<0) ) then RETURN([0,0,0]): fi:
S:=PROPSOL(A,B,C):
#note that  S[3] <> 0 unless  d1  is a square -- never true for us
LRT(5,`Done!`):
if FATAL then RETURN([0,0,0]): else
RETURN([bc*S[1],S[2],d0b*sqrt(d0a)*S[3]/sqrt(d0d)]); fi:
#need _not_ be integers, or in lowest terms!
end:

######################## NUMBER THEORY BASICS ###############################
#Mostly copied from Maple with interrupts to allow only limited factorization 

PROPSOL:=proc(a1,b1,c1) #used only in  NORM0, so presumably never needed
local sg,alpha,beta,gamma2,l,a,b,c,A,B,C,index,s,r,Q,ps;
#small modif of isolve/pyth/propsol,`Copyright 1993 by Waterloo Maple Software`
LRT(10,`Entering propsol`,args);
    l := sort([a1,b1,c1],proc(x,y) evalb(abs(x) < abs(y)) end);
    sg := signum(l[3]);
    a := l[1]*sg;
    b := l[2]*sg;
    c := l[3]*sg;
    index := abs(a*c);
    if index = 1 then
        if a+c = 0 then ps := [1,0,1]
        elif a+b = 0 then ps := [1,1,0]
        elif b+c = 0 then ps := [0,1,1]
        fi:
    else
        r := MSQRT(-b/a,c);if FATAL then RETURN([0,0,0]):fi:
        if abs(1/2*c) < r then r := c-r fi;
        Q := (a*r^2+b)/c;
        if Q = 0 then ps := [1,1,0]
        else
            A := igcd(b,c*Q);
            B := a*b/A;
            alpha := r/A;
            beta := b/A;
            C := Q/A;
            gamma2 := SQRPART(C);
            C := C/gamma2;
            s := procname(A,B,C);if FATAL then RETURN([0,0,0]):fi:
            ps :=
                [A*alpha*s[1]-beta*s[2],s[1]+a*alpha*s[2],C*sqrt(gamma2)*s[3]];
        fi;
    fi;
    ps := map(proc(x,y) abs(x/y) end,ps,igcd(op(ps)));
    if a = a1*sg then if b = b1*sg then ps else [ps[1],ps[3],ps[2]] fi;
    elif b = a1*sg then
        if a = b1*sg then [ps[2],ps[1],ps[3]] else [ps[2],ps[3],ps[1]] fi;
    else if a = b1*sg then [ps[3],ps[1],ps[2]] else [ps[3],ps[2],ps[1]] fi;
    fi:
end:

SQRPART:=proc(n) #used in NORM, PROPSOL, NORMXa, NORM0; uses IFACTOR
local i,f,m;   
#small modif of isolve/pyth/sqrprt  `Copyright 1993 by Waterloo Maple Software`
# NOTE: this may be just _part_ of the square part;  I think that any square
#  factors remaining will just slow down  propsol, but still give correct ans.
LRT(10,`Entering sqrpart`,args);
    m := 1;
    f := IFACTORS(n);
    for i in op(2,f) do  m := m*i[1]^(i[2]-(i[2] mod 2)) od;
    m
end:

MSQRT:=proc(x,p1) #used in PROPSOL NORM (NORMYa NORMXb NORMZb too); uses FACINT
local p,fact,y,res,bas,f,FFF,j;
global FATAL:
options remember;
#small modif of numtheory/msqrt `Copyright 1993 by Waterloo Maple Software`
LRT(10,`Entering msqrt`,args);
p:=abs(p1):
y := modp(x,p);
if y < 2 then RETURN(y) fi;
if p < 80 then
    for j to iquo(p,2) do  if modp(j^2,p) = y then RETURN(j) fi od;
    FATAL:=true:
    #this reallly shouldn't happen -- I only call MSQRT when sqrt exists!
    RETURN(FAIL)
fi;
if isprime(p) then RETURN(_msqrtprime(y,p)) fi;
res := 0;
bas := 1;
fact := FACINT(p); #I assume this has already tried to find every avail. fact.
CHKSNAG():if FATAL then RETURN(0): fi:
#Currently not really treated as "fatal", but it's "locally fatal"
for f in fact do
    for j while modp(p,f^(j+1)) = 0 do   od;
    _msqrtpk(modp(y,f^j),f,j);
    if " = FAIL then RETURN(FAIL) fi;
    res := numtheory[mcombine](bas,res,f^j,");
    bas := bas*f^j
od;
RETURN(res):
end:

################ QUARTIC REDUCTION ######################################
#New. I should meld with local solvability.
#Reduction as first published by B,S-D (I); verbatim in Cremona, Serf, etc.

INVARS:=proc(a,b,c,d,e): #computes  I  and  J. Note: Maple dislikes "I".
RETURN([12*a*e-3*b*d+c^2,72*a*c*e+9*b*d*c-27*a*d^2-27*e*b^2-2*c^3]):end:

IJREDUCE:=proc(a,b,c,d,e)
local s,t,i,a1,b1,c1,d1,e1,II,J,ff,g,ffold,g0,fp,ep,oo,qt:
global QTIME:
LRT(10,`Entering ijreduce`):
qt:=time():
g:=igcd(a,b,c,d,e):
s:=IFACTORS(g)[2]:
t:=1:for i in s do t:=t*i[1]^(i[2]-modp(i[2],2)):od:
a1:=a/t:b1:=b/t:c1:=c/t:d1:=d/t:e1:=e/t:g:=g/t:
ff:=a1*u^4+b1*u^3+c1*u^2+d1*u+e1:
oo:=[]:for i from 0 to 4 do oo:=[coeff(ff,u,i),op(oo)]:od:
INVARS(a1,b1,c1,d1,e1): II:="[1]:J:=""[2]:
s:=IFACTORS(II)[2]:
for i in s do 
 if (i[1]>3 and i[2]>3 and modp(J,i[1]^6)=0) or 
  (i[1]=3 and i[2]>4 and modp(J,3^9)=0) or
  (i[1]=3 and i[2]=4 and modp(J,3^6)=0 and modp(J,3^7)>0 and modp(4*II^3-J^2,3^15)=0) or
  (i[1]=2 and i[2]>5 and modp(J,2^9)=0 and modp(8*II+J,2^10)=0) then 
	ffold:=ff: #just in case I've messed this up
	if modp(ff,i[1])=0 then g0:=i[1] else g0:=1: fi:
	fp:=Factor(ff/g0) mod i[1]: #seems to be const or const*lin^4
#Nope! Try mod 3:
#                 4            3             2
#        - 47175 x  + 2317140 x  + 36635760 x  - 450415320 x - 683262720
	if degree(fp,u)>2 then #we assume it has at least a triple root
		fp:=factors(fp)[2]:
		if fp[1][2]>2 then fp:=fp[1][1] else fp:=fp[2][1]: fi:
		ep:=solve(fp,u):
		ff:=expand(subs(u=i[1]*u+ep,ff/i[1]^4)):
	else #assuming there's a triple root at infinity
		ff:=subs(u=u/i[1],ff):
	fi:
	oo:=[]:for i from 0 to 4 do oo:=[coeff(ff,u,i),op(oo)]:od:
	for i from 1 to 5 do if not type(oo[i],integer) then 
		WARN(`IJREDUCE got carried away; self-corrected.`);
		ff:=ffold:fi:od:
fi:od:
oo:=BSDREDUCE(op(oo)):
if WANTSEARCH then ADDREQUEST(op(oo)):fi:
QTIME:=QTIME+time()-qt:
RETURN(op(oo)):
end:

BSDREDUCE:=proc(a,b,c,d,e)
local f,numroots,ss,i,h,BSDp,BSDx,p,q,p2,q2,flag:
LRT(10,`Entering bsdreduce`):
Digits:=2*(length(a)+length(b)+length(c)+length(d)+length(e)+5):
f:=a*x^4+b*x^3+c*x^2+d*x+e:
numroots:=nops(realroot(f)):LRT(6,`There are`,numroots,`real roots`);
ss:=[evalf(allvalues(RootOf(f)))]:
flag:=+1:  # in one case we effect a sign change
#
#Create covariant quadratic. (Sure would be nice to do it exactly (Digits big!)
if numroots=2 then
	h:=1:for i in ss do if not type(i,realcons) then h:=h*(x-i):fi:od:
	h:=expand(h): 
elif numroots=0 then
	p:=-2*Re(ss[1]):
	q:=Re(ss[1])^2+Im(ss[1])^2;
	p2:=b/a-p:
	q2:=(e/a)/q:
	h:=expand((x^2+p*x+q)/sqrt(4*q-p^2)+(x^2+p2*x+q2)/sqrt(4*q2-p2^2));
else #numroots=4
	ss:=sort(ss):
	if ss[3]-ss[1] < ss[4]-ss[2] then
	 p:=+ss[1]+ss[3]:q:=ss[1]*ss[3]:p2:=+ss[2]+ss[4]:q2:=ss[2]*ss[4]:
	 flag:=-1: #note that we need to replace x by -x here.
	else
	 p:=-ss[2]-ss[4]:q:=ss[2]*ss[4]:p2:=-ss[1]-ss[3]:q2:=ss[1]*ss[3]:
	 flag:=+1:
	fi:
	h:=(p2-p)*x^2 + 2*(q2-q)*x + (p*q2-p2*q):
fi:
#
SHRINKQF(coeff(h,x,2),coeff(h,x,1)/2,coeff(h,x,0)):
BSDp:=expand(numer(subs(x=flag*("[4]*x+"[5])/("[6]*x+"[7]),f))):
BSDx:=[flag*(""[4]*x+""[5])/(""[6]*x+""[7])];
LRT(5,`The transformation used is x |-> `,op(BSDx));
LRT(5,sort(BSDp));
RETURN([coeff(BSDp,x,4),coeff(BSDp,x,3),coeff(BSDp,x,2),coeff(BSDp,x,1),coeff(BSDp,x,0)]);
end:

SHRINKQF:=proc(a,b,c) 
local dis,x,r,s,u,t,h,al,be,ga,de,blah,duh,v1,v2,a1,b1,c1:
#return a nice quad form integrally equiv to ax^2+2bx+c
#and a change-of-variables formula to get there.
LRT(10,`Entering shrinkqf`,a,b,c):
Digits:=10+2*(length(a)+length(b)+length(c)):
dis:=b^2-a*c:
	#note: we are allowing _improper_ equivalence too
	#(GL(2Z) not SL2Z, ie x|-> -x is assumed OK)
	#Also note: real coeffs OK here. Also: it's early in Cohen's book
if a<0 then a1:=-a:b1:=-b:c1:=-c:else a1:=a:b1:=b:c1:=c:fi:
duh:=evalf(sqrt(-dis));
blah:=evalf(sqrt(a1)):
lattice([[blah,0],[b1/blah,duh/blah]]):
v1:="[1]: v2:=""[2]:
be:=round(v1[2]*blah/duh):al:=round(v1[1]/blah-be*b1/a1):#v1=al*1st+be*2d
de:=round(v2[2]*blah/duh):ga:=round(v2[1]/blah-de*b1/a1):#v2=ga*1st+de*2d
h:=numer(subs(x=(al*x+ga)/(be*x+de),a*x^2+2*b*x+c)):
RETURN(coeff(h,x,2),coeff(h,x,1)/2,coeff(h,x,0),al,ga,be,de,`(pos) def`):
end:

#following not used, nor complete.

EQUIV:=proc(a,b,c,d,e,a1,b1,c1,d1,e1) #incomplete test for equivalence
local cr,qr,al,be,ga,de:
LRT(10,`Entering equiv`):
cr:=[evalf(allvalues(RootOf(a*x^4+b*x^3+c*x^2+d*x+e)))];
qr:=[evalf(allvalues(RootOf(a1*x^4+b1*x^3+c1*x^2+d1*x+e1)))];
al:=cr[1]*qr[1]*(qr[2]-qr[3])+cr[2]*qr[2]*(qr[3]-qr[1])+cr[3]*qr[3]*(qr[1]-qr[2]);
be:=cr[1]*qr[1]*(cr[2]*qr[3]-qr[2]*cr[3])+cr[2]*qr[2]*(cr[3]*qr[1]-qr[3]*cr[1])+cr[3]*qr[3]*(cr[1]*qr[2]-qr[1]*cr[2]);
ga:=(cr[2]*qr[3]-qr[2]*cr[3])+(cr[3]*qr[1]-qr[3]*cr[1])+(cr[1]*qr[2]-qr[1]*cr[2]);
de:=qr[1]*cr[1]*(cr[2]-cr[3])+qr[2]*cr[2]*(cr[3]-cr[1])+qr[3]*cr[3]*(cr[1]-cr[2]);
be:=be/al:ga:=ga/al:de:=de/al:al:=1:
RETURN([al,be,ga,de]);
end:

############## LOCAL SOLVABILITY OF QUARTICS ##############################
# Fast local solvability (after Siksek)
# For p=3,5 (,...?) it may be faster just to do  exhaustive search, as for p=2.
# Note: for _groups_ of quartics I should probably keep track of group 
# structure for each relevant   p. Not yet done.

LOCALLYSOLV:=proc(a,b,c,d,e)
local i,q,x,f,p,n,ls,bads: 
global LSTIME,LIKELIES,SNAG:
#
# Checks to see if  y^2=quartic  is everywhere locally solvable.
#
LRT(10,`Entering locallysolv`,args);
ls:=time():
SNAG:=0;
f:=a*x^4+b*x^3+c*x^2+d*x+e: #assumed integral,
for i from 0 to 4 do q[i]:=coeff(f,x,i):od:
#Warning: no checks now for infinite loops (should not take more than
#ordp(discrim(f,x),p) trips through  ZPSOLV. (maybe a little worse for p=2?)
if nargs=5 then LIKELIES:=[2]:else #else we've supplied the bad primes
	bads:=args[6] minus {-1}:fi:
for p in LIKELIES do
	#NOTE: it seems a few primes keep disallowing _many_ d3's for each d1.
	#      Obviously should check them first.
	LRT(6,`   ...checking quartic at p=`,p);
	if not QPSOLV(q[4],q[3],q[2],q[1],q[0],p) then 
	    PRT(5,`quartic is not solvable at`,p):
	    LSTIME:=LSTIME+time()-ls:
	    RETURN(false):fi:
	od:
if not REALSOLV(q[4],q[3],q[2],q[1],q[0]) then 
	PRT(6,`quartic is not solvable in the reals.`):
	LSTIME:=LSTIME+time()-ls:
	RETURN(false):fi:
if nargs=5 then bads:=FACINT(abs(discrim(f,x))):fi: #may set SNAG
for p in bads do if not member(p,LIKELIES) then
	LRT(6,`   ...checking quartic at p=`,p);
	if not QPSOLV(q[4],q[3],q[2],q[1],q[0],p) then 
	    LIKELIES:=[op(LIKELIES),p]:
	    PRT(5,`quartic is not solvable at`,p):
	    LSTIME:=LSTIME+time()-ls:
	    RETURN(false):fi:
fi:od:
#Well, I guess this quartic passes muster!
LSTIME:=LSTIME+time()-ls:
RETURN(true);
end:

REALSOLV:=proc(a,b,c,d,e)
local x:
LRT(10,`Entering realsolv`,args);
if a>0 then RETURN(true):fi:
if e>0 then RETURN(true):fi:
if nops(realroot(a*x^4+b*x^3+c*x^2+d*x+e))>0 then RETURN(true):fi:
#Hey, OK, so I'm too lazy to program sturm sequences for a quartic myself...
RETURN(false):
end:

TWOSOLV:=proc(a,b,c,d,e)
local x,i,j,f,f2,q:
global TROUBLE:#Something is wrong. I need to rewrite this. Try 1,0,-3,0,2 !
# Solvability at 2. If I goofed here, use the one in Cremona's book, it's fine
LRT(10,`Entering twosolv`,args);
TROUBLE:=TROUBLE-1:
if TROUBLE<0 then ERR(`TWOSOLV screwed up`):RETURN(true):fi:
f:=a*x^4+b*x^3+c*x^2+d*x+e;
for i from 0 to 7 do
	if modp(a*i^4+b*i^3+c*i^2+d*i+e,8)=1 then RETURN(true):fi: od:
#the only other way to solve y^2=f is to have f=0 mod 4. 
for i from 0 to 3 do if modp(a*i^4+b*i^3+c*i^2+d*i+e,4)=0 then
	f2:=expand(subs(x=i+4*y,f)):
	if f2=0 then RETURN(true):fi:#should probably ring bells too
	while modp(content(f2),4)=0 do f2:=f2/4: od: 
	for j from 0 to 4 do q[j]:=coeff(f2,y,j):od:
	if TWOSOLV(q[4],q[3],q[2],q[1],q[0]) then RETURN(true):fi:
	fi:od:
RETURN(false):
end:

QPSOLV:=proc(a,b,c,d,e,p):
LRT(10,`Entering qpsolv`);
if ZPSOLV(a,b,c,d,e,p) then RETURN(true):fi:
if ZPSOLV(e,d,c,b,a,p) then RETURN(true):fi:
RETURN(false):
end:

ZPSOLV:=proc(a,b,c,d,e,p)
local x,i,f,f0,f1,ff,sq,ep1,ep2,g,h,f_1,f_2,co:
global TROUBLE:#see TWOSOLV
#
#Quick method to test a quartic for solvability in p-adic integers.
#
LRT(10,`Entering zpsolv`);
if p=2 then TROUBLE:=100:RETURN(TWOSOLV(a,b,c,d,e)):fi:
f:=a*x^4+b*x^3+c*x^2+d*x+e:
f0:=(Factor(f) mod p):
if f0=0 then RETURN(STEP2(a,b,c,d,e,p)):fi:
if degree(f0,x)=0 then if L(f0,p)=1 then RETURN(true):else RETURN(false):fi:fi:
f1:=factors(f0):
ff:=f1[2]:sq:=true:
for i from 1 to nops(ff) do if (ff[i][2] mod 2)=1 then sq:=false:fi:od:
if sq=false then RETURN(true):fi: #so now assume it's const*square
if L(f1[1],p)=1 then RETURN(true):fi: #so now it's (nonsq const)*square
#the square factor is (irr quadratic)^2 or (lin)^4 or (lin)^2(lin^2) or lin^2
if nops(ff)=1 and degree(ff[1][1],x)=2 then RETURN(false):fi:#so it factors...
ep1:=solve(ff[1][1],x):                                  #...and here is...
if nops(ff)>1 then ep2:=solve(ff[2][1],x):               #...(or are) the roots
else ep2:=ep1:
 #Note: these primes are special: if searching for points on such a hgs sp.,
 #must have  x = ep1 mod p. So substitute x=ep1+p*x': quartic is "simpler".
 #Probably if return is "true" we should return  ep1. Likewise below
 #if precisely one ep_i yields a positive return from  ZPSOLV.
fi:
g:=(x-ep1)*(x-ep2): h:=expand((f-f1[1]*g^2)/p):
#if (subs(x=ep1,h)*subs(x=ep2,h) mod p) > 0 then RETURN(false):fi:
if (subs(x=ep1,h) mod p)=0 then
 f_1:=expand(subs(x=p*x+ep1,f)/p^2):
 for i from 0 to 4 do co[i]:=coeff(f_1,x,i):od:
 if ZPSOLV(co[4],co[3],co[2],co[1],co[0],p) then RETURN(true):fi:
fi:
if ep1=ep2 then RETURN(false):fi:
if (subs(x=ep1,h) mod p)>0 then RETURN(false): fi:
 f_2:=expand(subs(x=p*x+ep2,f)/p^2):
 for i from 0 to 4 do co[i]:=coeff(f_2,x,i):od:
 if ZPSOLV(co[4],co[3],co[2],co[1],co[0],p) then RETURN(true):fi:
RETURN(false):
end:
#
STEP2:=proc(a,b,c,d,e,p)
local x,i,f,f1,f2,f3,hasroot,nonrept,ep1,ep2,co:
LRT(10,`Entering step2`);
f:=a*x^4+b*x^3+c*x^2+d*x+e:
f1:=f/p:
if (f1 mod p)=0 then RETURN(ZPSOLV(a/p^2,b/p^2,c/p^2,d/p^2,e/p^2,p)):fi:
f2:=Factor(f1) mod p:f3:=factors(f2)[2]:
hasroot:=false:nonrept:=false:
for i from 1 to nops(f3) do 
	if degree(f3[i][1],x)=1 then 
		hasroot:=true:
		if f3[i][2]=1 then nonrept:=true:fi:
	fi: od:
if hasroot=false then RETURN(false):fi:  #so it's got a root...
if nonrept=true then RETURN(true):fi:    #...but all are repeated:f1=const*()^2
ep1:=solve(f3[1][1],x):ep2:=ep1:
if nops(f3)>1 then ep2:=solve(f3[2][1],x):fi:        #and these are the roots
for i from 0 to 4 do co[i]:=coeff(expand(subs(x=p*x+ep1,f1)/p),x,i):od:
if ZPSOLV(co[4],co[3],co[2],co[1],co[0],p) then RETURN(true):fi:
if ep1=ep2 then RETURN(false):fi:
for i from 0 to 4 do co[i]:=coeff(expand(subs(x=p*x+ep2,f1)/p),x,i):od:
RETURN(ZPSOLV(co[4],co[3],co[2],co[1],co[0],p)):
end:

################ NEXT-STEP DESCENT ? ####################################

TIDYQUART:=proc(a,b,c,d,e)
local f,g,h,b1,c1,d1,e1,co,i,su:
LRT(10,`Entering  tidyquart`,a,b,c,d,e);
f:=a*x^4+b*x^3+c*x^2+d*x+e:
g:=a*((x^2+b1*x+c1)^2+d1*(x+e1)^2):
h:=expand(f-g):
for i from 0 to 3 do co[i]:=coeff(h,x,i):od:
su:=solve({co[0],co[1],co[2],co[3]},{b1,c1,d1,e1})[1]:
g:=subs(su,d1):f:=sqrt(SQRPART(numer(g))):h:=sqrt(SQRPART(denom(g))):
LRT(6,a,`times the polynomial may be rewritten as`);
PRT(6,subs(su,(a*x^2+a*b1*x+a*c1)^2+(g*h^2/f^2)*(f*a/h*x+f*a/h*e1)^2));
LRT(6,`It factors over Q[ sqrt(`,-numer(g)/f^2*denom(g)/h^2,`) ].`);
#RETURN(subs(su,[a,b1,c1,d1,e1])):
RETURN(subs(su,[a,b1,c1,(g*h^2/f^2),(f*a/h), f*a/h*e1]));
end:

#The following lets us get more quartics equivalent (?) to the "other" d3 
#(hgs-sp +hgs-sp corresp to elemetn of order 2 in E'  ?) Interestingly inputs 
#which are PGL2Z-equiv need not give outputs which are PGL2Z-equiv?
#WARNING - BUGGY.....

REPEAT:=proc(a,b,c,d,e)
local AA,BB,TT,boo,ug:
LRT(10,`Entering repeat`);
TT:=TIDYQUART(a,b,c,d,e);
AA:=TT[1]: BB:=-TT[4]: 
#BB appears to be the value of "d1" attached to hgs sp. For T=19/2, d1=-1,
#the Other value is the same as (1-T^10)/5/d1 mod squares.
if not ((TT[4]=1) or (TT[4]=-30655331283885)) then 
	WARN(`Didn't get the expected splitting`);fi: #EXPERIMENTAL....
NORM(BB,AA);boo:=simplify([solve(X^2-BB=AA*("[3]/"[2]+m*(X-"[1]/"[2]))^2,X)]);
if type(boo[1],rational) then boo[2] else boo[1]:fi:
"-simplify((AA*x^2+AA*TT[2]*x+AA*TT[3])/(TT[5]*x+TT[6]));
numer(");discrim(",x);ug:="/SQRPART(content("));
IJREDUCE(coeff(ug,m,4),coeff(ug,m,3),coeff(ug,m,2),coeff(ug,m,1),coeff(ug,m,0));
RETURN(");
end:


FOLLOWUP:=proc(a0,b0,c0,d0,e0,d1)
local a,b,c,d,e,f,X,Y,Z,X0,Y0,Z0,L,R,nu,de,eq,dset,cset,lam,pr,p:
LRT(10,`Entering followup`,a0,b0,c0,d0,e0,d1):
X:=TIDYQUART(a0,b0,c0,d0,e0):
a:=X[1]:b:=X[2]:c:=X[3]:d:=X[4]:e:=X[5]:f:=X[6]:
#so we have a quartic a*y^2=(a*x^2+a*b*x+a*c)^2+d*(e*x+f)^2
R:=simplify((a*x^2+a*b*x+a*c)/(e*x+f)):#so R^2+d=a(y~)^2
X:=NORM(-d,a):X0:=X[1]:Y0:=X[2]:Z0:=X[3]:
#so get [X,Y,Z] with X^2+d*Y^2=a*Z^2. Then gen'l sol'n has X=
L:=simplify((X0*a*m^2-2*a*Z0*m+X0)/(a*m^2-1)/Y0):
#So we want L=R for some  x,m
dset:=FACINT(resultant(numer(L),denom(L),m)*resultant(numer(R),denom(R),x)):
#candidates for  d5. Equations to be solved are (the homog. forms of)
#numer(L)=d5*numer(R),  denom(L)=d5*denom(R).
#Combine into one over Q[S]. (Let's assume the hgs-sp "d1" is -1 for now...)
if not d1=-1 then WARN(`Procedures assume a (-1)-homogeneous space!`);fi: 
lam:=(1/e^2)*(b*a*e-2*f*a+2*a*S*sqrt((-e*b*f+f^2+e^2*c)/d1));
#S=alias(sqrt(d1))
if not type(sqrt((-e*b*f+f^2+e^2*c)/d1), rational) 
	then ERR(`Wrong splitting field?`):fi:
L:=numer(L)-lam*denom(L):X:=denom(icontent(expand(L)));L:=L*X:
#If I've chosen  lam  right, this rhs should be square*const..
#Compare coeffs of 1, S in sqrt(rhs/const) to find  x, y. Main action on LHS.
cset:=ODDINT(X*coeff(numer(R),x,2)) minus dset:
pr:=1:for p in cset do pr:=pr*p:od:
eq:=expand(simplify(subs(m=m/n,L)*n^2))=pr*d5*W^2:
X:=coeff(lhs(eq),m,2):Y:=coeff(coeff(lhs(eq),m,1),n,1):Z:=coeff(lhs(eq),n,2):
#eqn is  (2Xm+Yn)^2+(4XZ-Y^2)n^2=(4X*pr*d5)W^2
for L in powerset(dset) do
	d:=1:for p in L do d:=d*p:od:
	R:=POLAR(1,subs(S=I,4*X*Z-Y^2),subs(S=I,-4*X*pr*d));
	if GAUSSHASS(op(")) then print(R):fi:
od:
print(`Here is the single quadratic to be solved over Q[ sqrt(`,d1,`) ]`);
print(eq);
print(`d5's in`,dset);
end:

POLAR:=proc(A,B,C) #reduce a conic over Q[i]
local A1,B1,C1,A2,B2,C2,X:
LRT(10,`Entering polar`,A,B,C):
X:=GAUSSGCD(A,GAUSSGCD(B,C)):
A1:=A/X:B1:=B/X:C1:=C/X:
A2:=GAUSSGCD(B1,C1):
B2:=GAUSSGCD(C1,A1):
C2:=GAUSSGCD(A1,B1):
RETURN([A1/B2/C2*A2,B1/C2/A2*B2,C1/A2/B2*C2]);
end:

GAUSSHASS:=proc(A,B,C) #Tests for local solvability.
if not MINITEST(A,B,C) then RETURN(false):fi:
if not MINITEST(B,C,A) then RETURN(false):fi:
if not MINITEST(C,A,B) then RETURN(false):fi:
RETURN(true):
end:

MINITEST:=proc(A,B,C)
local i,p,fa,RE,IM:
fa:=GAUSSFACT(A)[2]:
for i in fa do
	if Im(i[1])<>0 then
		RE:=Re(i[1]);IM:=Im(i[1]):p:=RE^2+IM^2:
		subs(I=-RE/IM,[B,C]) mod p:
		if L(-B*C,p)=-1 then RETURN(false):fi:
	else
		-B*C mod i[1]: #needs to be a square.
		#x=y^2 iff 1= x^(p^2-1)/2 =x^(1+p)*(p-1)/2, so iff N(x)=[]
		CNORM("):
		if L(",i[1])=-1 then RETURN(false):fi:
	fi:
od:
RETURN(true):
end:

#discuss local solv? quartic local solv? find param'n?

#################### GAUSSIAN ARITHMETIC ################################
# Sort of temporary - I really need arith, over general quartic fields. #

CNORM:=proc(x):RETURN(expand(x*subs(I=-I,x))):end:

CPRIMEDEC:=proc(x)
local R, IM:
LRT(10,`Entering cprimedec`,x):
if modp(x,4)=3 then RETURN(x):fi:
NORM(-1,x):R:=abs("[1]):IM:=abs(""[2]):
if R>IM then RETURN(R+I*IM) else RETURN(IM+I*R) fi:
end:

GAUSSGCD:=proc(n0,n1) #GCD calculation in gaussian integers.
local x0,y0,x1,y1,g:
LRT(10,`Entering gaussgcd`,n0,n1):
x0:=Re(n0):y0:=Im(n0):x1:=Re(n1):y1:=Im(n1):g:=1:
while g<>0 do
 g:=DIVALGa(x0,y0,x1,y1,-1);
 x0:=x1:y0:=y1:x1:=g[1]:y1:=g[2]:
 g:=x1+I*y1:
od;
RETURN(x0+I*y0):
end:

GAUSSFACT:=proc(x) #factor the _gaussian integer_  x.
local ans,z,n,fa,i,j,k,p,pbar,pos,neg:
LRT(10,`Entering gaussfact`,x):
ans:=[]:
z:=x:n:=CNORM(x):
fa:=IFACTORS(n)[2]:
for i in fa do
  if modp(i[1],4)=3 then 
	k:=i[2]/2 : ans:=[op(ans),[i[1],k]]: z:=z/i[1]^k:
  else 
	p:=CPRIMEDEC(i[1]):pbar:=subs(I=-I,p):pos:=0:neg:=0:
	for j from 1 to i[2] do
		if denom(icontent(z/p))=1 then 
			z:=z/p:pos:=pos+1:
		else z:=z/pbar:neg:=neg+1: fi:
	od:
	if pos>0 then ans:=[op(ans),[p,pos]]:fi:
	if neg>0 then ans:=[op(ans),[pbar,neg]]:fi:
  fi:
od:
RETURN([z,ans]);#first component should be one of the 4 units.
end:

#################### INTEGER FACTORING ROUTINES ########################
# The only part of the whole program which isn't provably polynomial-time
#is the factorization of a numbers. Current code does very little factoring,
#but these routines are retained.
# We keep a list  KNOWNPR  of every relevant prime, as well as another
#list  UNKNOWN  of composite numbers, each prime to all of
#KNOWNPR as well as to each other; each is also not a square or cube, and is
#not divisible by any  p<1700. We maintain these attributes as the lists grow.

FACINT:=proc(n)
local z, S, i:  #        n   is assumed integral!
LRT(10,`Entering facint`,n);
z:=IFACTORS(n): #note - don't chksnag -- leave that to calling proc if needed.
S:={}:
if z[1]=-1 then S:={-1}:fi:
for i from 1 to nops(z[2]) do S:=S union {z[2][i][1]}: od:
RETURN(S):
end:

ODDINT:=proc(n)
local z, S, i:
LRT(10,`Entering oddint`,n);
z:=IFACTORS(n):
S:={}:
if z[1]=-1 then S:={-1}:fi:
for i from 1 to nops(z[2]) do if (z[2][i][2] mod 2) = 1 then S:=S union {z[2][i][1]}: fi: od:
RETURN(S):
end:

QUICKFACT:=proc(n)
local z,p,L,k,NEW_UNK,NOT_UNK,S,i:#check out factors already known to this run.
global UNKNOWN, KNOWNPR:
LRT(10,`Entering quickfact`,n);
z:=n:
L:=[]:
for p in KNOWNPR do
	k:=0:
	while (z mod p)=0 do z:=z/p:k:=k+1:od: 
	if k>0 then L:=[op(L),[p,k]]; fi:
	od:
NOT_UNK:={}:NEW_UNK:={}:
for p in UNKNOWN do
	k:=gcd(p,z): if (k>1) and (k<p) then 
		LRT(5,`Hah! factored an old composite a little:`,p):
		NOT_UNK:=NOT_UNK union {p}:
		S:=BREAKUP(p,z):
		for i in S do 
		 if isprime(i[1]) then KNOWNPR:=KNOWNPR union {i[1]}:
		   if i[3]>0 then L:=[op(L),[i[1],i[3]]]: z:=z/i[1]^i[3]:fi:
		  else NEW_UNK:=NEW_UNK union {i[1]}: fi: 
		od:
	fi:od:
UNKNOWN:=(UNKNOWN minus NOT_UNK) union NEW_UNK:
RETURN([L,z]): end:

BREAKUP:=proc(a,b)
local k,S,T,s,t,expa,expb,expk,basea,baseb,basek,ans:
#
# Given two integers (a,b) return a set of triples [ui,ei,fi] with
# the ui rel. prime and  a = Prod(ui^ei), b=Prod(ui^fi).
#
LRT(10,`Entering breakup`,args);
if a=1 and b=1 then RETURN({}):fi:
if a=1 then RETURN({[b,0,1]}):fi:
if b=1 then RETURN({[a,1,0]}):fi:
k:=gcd(a,b):
if k=1 then RETURN({[a,1,0],[b,0,1]}):fi:
ans:={}:
ROOTS(k):basek:="[1]:expk:=""[2]:
ROOTS(a/k):basea:="[1]:expa:=""[2]:
ROOTS(b/k):baseb:="[1]:expb:=""[2]:
#so you see how to recoup a,b. Note gcd(basea,baseb)=1
S:=BREAKUP(basea,basek):
for s in S do 
 if s[2]>0 then ans:=ans union { [s[1],expa*s[2]+expk*s[3],expk*s[3]] }:
 else T:=BREAKUP(s[1],baseb):
  for t in T do 
   ans:=ans union {[t[1],expk*t[2]*s[3],expb*t[3]+expk*t[2]*s[3]]}:
   od:fi:od:
RETURN(ans):
end:

ROOTS:=proc(n)
local x,m:
# Many unfactored numbers are perfect powers...
LRT(10,`Entering roots`);
if abs(n)=1 then RETURN([n,1]):fi:
if issqr(n) then x:=ROOTS(isqrt(n)):RETURN([x[1],2*x[2]]):fi:
m:=iroot(n,3):if n=m^3 then x:=ROOTS(m):RETURN([x[1],3*x[2]]):fi:
#repeat last line for 5th powers etc? -- Nah, not for elliptic curves!
RETURN([n,1]);
end:

IFACTORS:=proc(n)
local sg,nn,z,a,b,i,L,M,S,ft,ex:
global UNKNOWN,KNOWNPR,SNAG,FACTTIME:
#
# Plug-in replacement for  ifactors  but looks for current good primes,
# and has the good sense to quit after a while.
#
LRT(10,`Entering ifactors`,n);
if n=0 then RETURN([0,[]]);fi:
ft:=time():
SNAG:=0:sg:=sign(n):z:=1:
nn:=abs(n):L:=[]: 
if not factorproc=`external` or nn<FACTLIM then
 ifactor(nn,easy): #First handle all the small factors
 for i in convert(expand("),set) do if type(i,integer) then z:=z*i:fi:od:
 z:=gcd(nn,z):# try n:=(nextprime(10^17)*nextprime(10^18))^2 to see why needed
 L:=ifactors(z)[2]:# a workaround for the nonexistent  "L:=ifactors(z,easy):"
 nn:=nn/z:  #Now  nn  has no small factors. Any known factors?
 fi:
if nn>1 then
 QUICKFACT(nn):
 a:="[1]:nn:=""[2]:L:=[op(L),op(a)]:
 fi: #   nn has no small or known factors. Is it a power?
if nn>1 then
 ROOTS(nn):nn:="[1]:ex:=""[2]:
 else ex:=1:fi: #  nn  has no small or known factors nor roots. Prime?
if isprime(nn) then
 L:=[op(L),[nn,ex]]:nn:=1:fi: # Running out of tricks. Get fancy?
if nn>1 then
 if not factorproc=`external` then
	LRT(5,`Please wait while I factor`,nn):
	if factorproc=`maple` then a:=ifactors(nn)[2]
		else z:=ifactor(nn,factorproc):a:=ifactors(z)[2]:fi:
	LRT(5,`Thank you for waiting.`):
 else #ie if factorproc=`external`  
  if nn<FACTLIM then  #      I'll only wait a minute or so!
	LRT(5,`Using maple to factor a moderate bignum:`,nn):
	a:=ifactors(nn)[2]:
	LRT(5,`DONE!`);
  else 
	#no small factors, no old factors, no roots, not prime, not small, 
	#	not `maple`: GIVE UP!
	LRT(5,`I cannot factor(ize)`,nn);
	UNKNOWN:=UNKNOWN union {nn}:
	a:=[]:SNAG:=nn:
	fi: #end of: whether big or not
 fi: #end of: whether external or not
 L:=[op(L),seq([a[i][1], a[i][2]*ex], i=1..nops(a))]:
fi: #end of: do we need to get fancy?
b:={}: for i from 1 to nops(L) do (L[i]): b:=b union {L[i][1]}: od:
KNOWNPR:=KNOWNPR union b:
FACTTIME:=FACTTIME+(time()-ft):
RETURN([sg,L]):  #this is _NOT_ a complete set of factors unless SNAG=0
end:

CHKSNAG:=proc()
global SNAG,FATAL,FATALSET:
if SNAG>0 then ERR(`I cannot proceed without knowing factors of`,SNAG):
FATAL:=true:FATALSET:=[op(FATALSET),SNAG]:fi:
end:

############################ HELP FILES ################################

HELP:=proc():
lprint(`This program is designed to perform the second 2-descent on homogeneous`);
lprint(`spaces of the curves  y^2=x*(x^2+A*x+B). Type  MAXRANK(A,B)  and you`);
lprint(`should get an answer indicating the maximum rank the curve may have`);
lprint(`consistent with first and second 2-descents.`);
lprint(`A and B must be rational; if non-integral they will be converted.`);
lprint();
lprint(`Type-checking is performed only for the front-end routines MAXRANK,`);
lprint(`ONECURV, ONED1, and MANYD1. (Type e.g. ONED1();  for help.) Of course`);
lprint(`internal routines can be called by the user but should be used with care.`);
lprint();
lprint(`The program will explain something of what it is doing if you set`);
lprint(`     yaklevel:=4: `);
lprint(`(the default); setting  yaklevel  higher or lower will cause`);
lprint(`more (resp. less) information to be displayed.`):
lprint(`The source code is intended to be readable too (more or less).`);
lprint(`For a discussion of the theory involved visit`);
lprint(`            http://www.math.niu.edu/~rusin/research-math/descent/`);
lprint();
lprint(`I have added some efficient procedures for finding points on conics and`);
lprint(`for wading through large collections of possible conics and divisors`);
lprint(`Peak efficiency was not the goal, however, and clumsy tools were used`);
lprint(`whenever I spotted a way to get Maple to supply answers with a minimum`);
lprint(`of programming time (e.g.  factor (mod p), realroot, msqrt, jacobi,...)`);
lprint(`The principal exception was  ifactor, which caused early versions of the`);
lprint(`program to choke on the huge integers produced. This is no longer a problem.`);
lprint(`(Note that the program will look for and create a list prims.m of good primes)`);
lprint();
lprint(`If a homogeneous space "d1" passes the next descent test, a quartic will`);
lprint(`be displayed a solution of which induces a point on the original homogeneous`);
lprint(`space. No search routine is provided; Cremona's  ratpoint  works quickly`);
lprint(`for this purpose. An input file  ratpoint.in  will be created for it.`);
lprint(`Note that conceivably the homogeneous space "d1" may have points even though`);
lprint(`the printed quartic "d3" is not one of the ones which has solutions.`);
lprint();
lprint(`While I have tested the procedures used here, the usual disclaimers apply`);
lprint(`regarding bugs. Please send reports of difficulty to`);
lprint(`Dave Rusin, rusin@math.niu.edu`, dat);
lprint(`(Please include the above date in correspondence.)`);
lprint(`It is in particular true that for curves with three 2-torsion points,`);
lprint(`the rank reports are inaccurate. The routine  ONED1  is still OK.`);
end:

print(`Version of`,dat,`is ready.`);

#Recently noticed bugs:
#TWOSOLV can go to infinite loop (q.v.)
#FINDD1S adds useless divisors when there's 2-torsion.
#MAXRANK(0,-157^2) dismisses the useful hgs space! (Note: MWRANK works with
#	the altered curve c = 471, d = 49298).
#ONECURV messes up rank counts with TORSION=2 if e.g. y^2=x(x-1)(x-4); what
#	is ODDINT(4), eh?