This file will be a nice summary, I hope, of methods for solving "norm
equations"  x^2 - D y^2 = N  for rational  x, y,  as needed in the
iterated descent program. The problem is more classically stated as
requiring solutions to
                a x^2 + b y^2 + c z^2 = 0  
The issue here is whether a solution can be computed in polynomial time,
if the initial  a, b, and  c  are presented in factored form.

Until I rewrite this in an easily-read format, I'll just attach excerpts 
from two letters I wrote recently, summarizing what I knew.

==============================================================================

I have been writing a computer program to perform higher descents on
elliptic curves of the form  y^2 = x ( x^2 + A x + B ). One step in this
process is to find a rational point on the curves  y^2 = d1 x^2 + A x + d2
where  d1 d2 = B. Of course this is equivalent to an equation of 
the form (*). The connection is explained, for example, in a paper by
Birch and Swinnerton-Dyer (Jour. Reine Angew. Math. v218 1965). They 
state that solving equation (*) is precisely the difficulty in
computation which they experienced.

I attach below a summary of methods I know for studying equation (*).
In general, the procedures are insufficient if  a, b, c  are very large.
I have found some simple methods which appear to be efficient but I
have not proved that they always succeed in polynomial time. I would 
like to know if you have any ideas which would allow me to solve (*)
more rapidly. Can you also suggest other work which has been
done in recent years?

Cremona explained that you have been studying (*) over general number 
fields. I hope that means the solution over the rational field  Q  is easy!

      O /
=======X=====COUPEZ ICI=============================================
      O \

Here is what I know:

1. Equation (*) is related to factorizing. Since we do not know an efficient
method to factorize large integers, no efficient solution to (*) exists
in general. (If we could solve (*) without a factorization of  a, say,
then we could compute a square root of  (-b)  mod  a. With a large number
of square roots known mod  a, we would eventually find two numbers  x, y  with
x^2=y^2 mod  a  but   x+y <> 0 mod  a   and  x-y <> 0 mod  a.  Then
(x+y)  has a proper common factor with  a.)
	So if we assume that all large numbers are hard to factorize, then
equation (*) is hard to solve.

2. On the other hand, if we assume that all large numbers are easy to
factorize, then (*) is easy to solve. Of course we can use the Hasse principle
to determine whether (*) has a solution or not. In addition, we could
use the same procedure to _compute_ a solution to (*) _modulo  a*b*c_.
(Here I use the fact that square roots mod p are easy to compute.)
	One may then use the  LLL  lattice reduction procedure on a
sublattice of  R^3  to find a vector on which  a x^2 + b y^2 + c z^2  is a
_small_ multiple  k*a*b*c  of  a*b*c. If that multiple is zero, we are 
done. If this  k=-1, then there is a standard reduction which produces a
solution to  (*). But in general, I do not know a way to find a solution
to (*) starting only from a solution to (*) modulo  a*b*c.  Is there a
method known?
	(I do know that if (*) is solvable then there is a "small"
solution:  one with  |x|^2 < |b*c| etc. Thus the  k  above can be made
very small. But even if  k=1  I do not know how to produce a solution 
of (*) from a solution of
	a x^2 + b y^2 + c z^2 = a b c.  )

3. There is an alternative solution if factorizations are easy. I think
this procedure goes back to Gauss.  Equation  (*)  is equivalent to an
equation of the form 
	x^2 - D y^2 = N.
So the algorithm is easy: solve  x^2 - D y^2 = 0  mod  N;  this gives
x^2 - D y^2 = k N  for some  k < N. By induction we can solve
x^2 - D y^2 = k. Then we form a quotient in the field  Q[ D^(1/2) ]
and find an element with  norm = x^2 - D y^2 = N.
	Unfortunately, this requires solving  X^2=D  mod  k, which
requires factorizing  k.  This is very hard in practice if  k  is large.
	This is the procedure used in Maple in the procedure  isolve.


So  (*)  is hard if factorizations are all hard, and it is easy if
factorizations are all easy. However, my case is  intermediate: in the
original equation (*), the terms  a, b, c  _are_ large but their
factorization is known. Is it possible to solve (*) in this case
without needing to factor other large integers?

4. Here is the simplest example: is there a fast procedure to solve
	x^2 + y^2 = p
if  p  is prime (p=1 mod 4)? In this case we know there is an _integral_
solution. It turns out that there _is_ a fast solution, due to Rabin
(it is in a 1977? paper of Rabin and Shallit). If  x0^2 + y0^2 = 0 mod p,
then a prime divisor of  p  divides  x0 + y0 i in  Z[i]. Of course that
prime divides  p + 0 i  too, so we may use the division algorithm
to find a generator   x + y i  of the ideal  (p, x0 + y0 i). Then
x^2 + y^2 = p.
	This procedure may be used if  p  is not prime, too (if we can
factorize  p). We may try to generalize to the case  x^2 - D y^2 = N,
but  Q[ D^(1/2) ]  has no easy division algorithm, for almost all  D.
However, a simple generalization of the procedure seems to work in
many cases. That means many fewer numbers must be factored. Do you know
of any work which would explain why this method can be as successful 
as it is? (Of course for general  D  there are number-theoretic reasons
why  x^2 - D y^2 = N  could be solvable for _rational_  x  and  y  but
not for _integral_  x  and  y. That is  OK  for me.)

5. In the example above, we do not need  p  to be prime (if we know its
factors). There is a method described in Cohen's book to solve
x^2 - D y^2 = p  for general  p, but I think in this case we _do_
need  p  to be prime. Is there an easy extension of that procedure
to the case of composite  p  whose factors are known? It is certainly
possible that I have missed some easy generalization. (I cannot find my copy
of Cohen's book so I cannot remember the name of the procedure, sorry.)


Let me tell you what I have used so far. It seems to work quickly in all
the cases I have tried. It does not require factorization if the initial
coefficients  a, b, c  are already factorized. It does not produce the
solution [x:y:z]  of smallest height, but the solutions are fairly small.
However, I do not know that this will always be true. I do not even
know that the procedure will always succeed. Do you recognize it? Is
there a proof that it will work? Can it be made more efficient?

The trick is to use the factorization of  a,b,c  only to compute a
number  M  such that  
(**)	M^2 = -a b  mod  c,  M^2 = -b c mod a,  M^2 = -c a mod b
(Such an  M  is a "certifcate" that  (*)  is solvable). 
	Proceeding as in (4) above, we find a solution to
(***) a x^2 + b y^2 + k c z^2 = 0  for some small  k. Then it is sufficient to
solve a related equation which now has  k  as a coefficient. But we
do not need to factor  k.  Instead, we need only find another integer  M
which satisfies the corresponding condition (**). This is not hard
using both  M  and the solution to (***). Therefore we may proceed
inductively.
	(There are some difficulties since we cannot assume the 
corresponding  a,b,c  at each stage are square-free, since we do not
factorize them. I am finishing that problem now, I hope.)

==============================================================================
(I was advised to consult a paper by Pollard and Schnorr, IEEE, 1985?)
==============================================================================

Thank you for your prompt response to my letter. Your recommendation
to look at the paper of Pollard and Schnorr led to unexpected good ideas.

Although the paper was not directly of interest to me, reading it
helped me understand what I was already doing, and convinced me that
the LLL routine (which I mentioned in my previous mail) should be
perfect for this application. (A colleague, Neil Dummigan, had
suggested it to me weeks ago but it has taken me until now to believe
that this is the right approach.)

I don't know if this is really of interest to you but I'll summarize
below what I now understand. (It helps me to explain it; also I hope this 
is of direct use to John Cremona.) I apologize for the length of this letter.

You wrote:
>All this is probably not the answer you expected, but
>it is what I know on the subject.

Indeed, this was most helpful; thanks again for the reference. I would
appreciate a copy of your current paper when you finish it.

dave

======================================================================

As I noted in my last letter, my problem may be expressed as: find a 
solution to the equation (*) x^2-Dy^2=N  assuming a factorization of  N.

The paper by Pollard and Schnorr is interesting. They show how to solve
(*) x^2-Dy^2=N  in integers mod  M. Now, I wanted to solve this with
rational  x  and  y, and you prefer to study integer  x  and  y. Of course
it is clear that the three problems are different yet related; it is by
thinking about the differences that I have understood what is going on.

First let me comment that if we compare the Pollard and Schnorr solution
to others I mentioned last time, we see almost all of them follow these steps:

0) Use a simple identity, if necessary, to allow the assumption that
	|D| < |N|.
1) Obtain a solution  x^2 - D y^2 = k N  for some  k  (a solution which is
	nontrivial in some way, e.g., GCD(x,N)=1 )
2) Use (1) to obtain a solution of  x^2 - D y^2 = k N  for a much smaller  k
	(roughly k=sqrt(|D|), for example)
3) Find (by induction) a solution to  x^2 - k y^2 = D and combine with the 
	solution in step (2)  by some simple identities to obtain a 
	solution to (*).


Step (1) is about as hard to do in integers as factoring  N; for me this is
no problem as  N  (and D) are factored, but in step(3), k is not factored, so
I wish to avoid returning to step (1) at that time.  All the methods I had 
seen simply assume step (1) could be carried out easily (that is, more or 
less, by assuming  N  is factored). 

Pollard and Schnorr however have a new approach to solve step (1) mod  M:
just choose a prime  p  with  p = N mod M and with  D  a square mod  p; 
then it is easy to solve  x0^2 - D = k p  in integers and reduce mod  M.
Note that they do not need to factor  M, N, or D. They have another method 
which is provably efficient (assuming the Generalized Riemann Hypothesis): 
select random  x  and  y  mod  M until  p=( k*(x^2-Dy^2) )mod M  is prime.  


It is in step (2) that the methods I mentioned last time differ.
But now I see they can all be expressed in the same language.
Indeed, consider the lattice
	L = { (u,sqrt(|D|)*v):  u*y0 = v*x0  mod N }
in  R^2, where  (x0,y0) is the solution found in step (1). For all points
(u,sqrt(|D|)*v) in the lattice we have y0^2(u^2-Dv^2)=v^2(x0^2-Dy0^2)=0 mod N 
and so (assuming gcd(x0,N)=1)  u^2-Dv^2 is a multiple of  N. On the other hand,
the lattice must contain a point in the disk with area  4|N|sqrt(|D|),
for which  u^2+|D|v^2 < (4/pi sqrt(|D|)) |N|. We can view step (2) as asking
for a quick way to compute such a short vector in  L.

(2A) I guess today the most efficient is to use the LLL basis reduction
routine. With it we may find a vector in  L  with  u^2 + |D| v^2 at most
twice this size, so that  u^2 - D v^2 = k N  with  |k| < 8/pi sqrt(|D|).

(2B) Pollard and Schnorr's approach is essentially Gauss's cyclic reduction:
repeatedly replace a basis {v1,v2} of  L  with an optimal {v2, v1 + q*v2}.
At least that's how I read Disquisitiones: Gauss prefers to think of the
basis of  Z^2  as fixed and the quadratic form as changing, but that's
equivalent to expressing a fixed quadratic form on different bases of  L.

(2C) Cornacchia's algorithm [I found my copy of Cohen's book!] can be 
described as searching through  L  using the Euclidean algorithm only on
the first coordinates ("u") of the elements of  L. Since the focus of that
algorithm is to find integer solutions, if the final vector chosen in this
way does not have  u^2-Dv^2=N  then for prime  N  the algorithm halts;
there is no step(3). 

(2D) The division algorithm I mentioned last time can also be viewed in
lattice terms! Note that the lattice is closed under the endomorphism
T(a,sqrt(|D|)*b) = (D*b, sqrt(|D|)*a); the division algorithm in Z[sqrt(D)]
simultaneously selects the optimal integers  q  and  r  for replacing a basis
{v1,v2} with a basis {v2, v1 + q*v2 + r*T(v2) }. Thus in practice it
should run even more quickly than Gaussian reduction.

(2E) Maple's isolve is fastest in step(2): it does nothing! (Of course its 
"k" is not nearly as small as  sqrt(|D|) but at least we have  |k| < |N|.)


In step (3), all methods return to step (1). This is slow if  k  must
be factored, but I described a method at the end of my last letter
which avoids factorization. (Pollard and Schnorr avoid it too, as I noted 
in step(1) above). Thus there certainly exists a polynomial-time algorithm
to solve (*) rationally if  N  and  D  are factored.


Now, why do all the algorithms take step(2) only as far as  k=sqrt(D)?
It is because that is the best one can do when using _integer_  x  and  y,
as we have discussed in previous letters. So if you are accustomed to
solving integrally, step(2) is a natural middle step. But for my
purposes, this is inappropriate. I should write the problem homogeneously
(as I did when I first wrote to you) and find integer solutions to  
	(**)	x^2 - D y^2 = N z^2, 
treating  z  equally with  x  and  y. I see now that the LLL approach 
which I mentioned last time should do precisely that.


You mentioned that you didn't understand what I meant by this approach (again
maybe this is because it is useless if you want _integral_ solutions to (*).)
I'll explain it in the same terms as steps (0)-(3) above, but in fact this
is simply a common proof of the Hasse principle for quadratic forms
(see e.g. Cassels' Geometry of Numbers for technical details).

0') In the symmetric setting it is not necessary to exchange  D  and  N.

1') Instead of solving  x^2 - D y^2 = 0  mod  N, we suppose we can solve
equation (**) mod  D*N. Then we have square roots  a = sqrt(D) mod N  and  
b = sqrt(N) mod D. (Computing all this is easy in my setting). 

2') Recall that  step(2) asked for the shortest element in a Z^2 lattice.
This time let  L'  be the lattice
	L' = { (u,sqrt(|D|)*v, sqrt(|N|)*w):  u=v*a mod N, u=w*b mod D }
(There should also be 2-adic congruences, and the conditions should be
modified if  N  and  D  are not relatively prime.) Then each element of
L'  determines integers  u, v, and  w  with  u^2 - D v^2 - N w^2 = 0 mod  D*N. 
On the other hand, there is surely a nonzero element in  L'  with  
	u^2 + |D| v^2 + |N| w^2 <=  k*|DN| 
for some small constant  k. Thus for this short vector, (u^2 - D v^2 - N w^2)
will be a small multiple of  D N. We set  step(2') to be the task of finding 
such a vector.

Gauss's method does not apply to lattices in  R^3  but Vallee has an algorithm 
to compute the shortest vectors in  L' (see Math Rev. 88h:11097 - I haven't
seen her paper.) Now, if the lattice  L'  is properly defined with 2-adic
congruences, then the shortest vectors have  k=0  ("k<3/4") by a theorem of 
Holzer (reproved by Mordell), so Vallee's procedure would give another 
efficient method for solving (**) this time giving small solutions.
(see example below).

We can still use  LLL  to look for short vectors in  L'; at worst, the vector 
found by  LLL  will have length too large by a factor of  4, so for the
short vector it finds we should have  |k|=0,1,or 2. Even |k|=2 seems
impossible as I read Cohen pp 84-85 ("part(3)").

3') If  k=0, we are done. If  k=-1, there is an identity (see e.g.
Borevich and Shafarevich p.64) which constructs a solution to (**).
This would be necessary for e.g. x^2 - 13 y^2 = 101 z^2 which leads to k=-1.
But there does not seem to be any identity to complete step(3) if for 
example k=+1. This is why I was reluctant to pursue this technique before.

What seems to be true is that if  L'  is constructed with the 2-adic 
congruences too, we _always_ have  k=0  for the vector found by LLL 
(i.e. no step(3') is needed). I can't prove that yet; I'm hoping 
that |k|=1 can be excluded by a parity argument.


Let me illustrate the success of the  LLL  method with this example I 
encountered in my computations: solve
	x^2 -123452853 y^2 = 375556542871687561 z^2.

Maple's  isolve  looks for all integers solutions; a simple choice of 
parameters in the presented solution yields
 x = 1046604075188381761638228306282085675127460591395491714850954855412939829
 y = 94195864965456602377931994957413486821431595250953785393138248206832
 z = 485077328854304825807066797821006059543420627128191447388327

If instead I spend more time in step(2) I can find a minimal vector
in the 2-dimensional lattice  L', namely
	 [11498644506, 3840749]
for which  u^2 - D v^2 = -4497*N. (This takes about a dozen iterations of
Gaussian reduction or the Euclidean algorithm) If one repeats
steps (1)-(2) three more times one arrives at a rational point
 x = 953325956294504993
 y = 55568067874461
 z = 1185298274

If instead we ask  LLL to find a short vector as in step (2') it finds
 x = 3983425934511
 y = 335173812
 z = 2307
which happens to have  x^2-Dy^2=Nz^2, i.e., k=0  and no  step(3')  is needed.

Clearly this example shows it is preferable to use  LLL  at least if
it can find a point with   k=0. If not, it might be worthwhile to
encode  Vallee's procedure.


So you see that reading  Pollard and Schnorr  was wise not because I
wanted to solve  mod  M  but because it forced me to think about the
difference between solving rationally and solving integrally. I am
glad you recommended it.

dave