*** Reduction Condition:

  -- LLL: the classical LLL reduction condition.

  -- BKZ: Block Korkin-Zolotarev reduction.
     This is slower, but yields a higher-quality basis,
     i.e., one with shorter vectors.
     See the Schnorr-Euchner paper for a description of this.
     This basically generalizes the LLL reduction condition
     from blocks of size 2 to blocks of larger size.


************* Calling Syntax for LLL routines ***************

long 
[G_]LLL_{FP,QP,XD,RR} (mat_ZZ& B, [ mat_ZZ& U, ] double delta = 0.99, 
                       long deep = 0, LLLCheckFct check = 0, long verbose = 0);

* The [ ... ] notation indicates something optional,
  and the { ... } indicates something that is chosen from
  among several alternatives.

* The return value is the rank of B (but see below if check != 0).

* The optional prefix G_ indicates that Givens rotations are to be used;
  otherwise, classical Gramm-Schmidt is used.

* The choice FP, QP, XD, RR determines the precision used.

* If the optional parameter U is given, then U is computed
  as the transition matrix:

     U * old_B = new_B

* The optional argument "delta" is the reduction parameter, and may
  be set so that 0.50 <= delta < 1.  Setting it close to 1 yields
  shorter vectors, and also improves the stability, but increases the
  running time.  Recommended value: delta = 0.99.

* The optional parameter "deep" can be set to any positive integer,
  which allows "deep insertions" of row k into row i, provided i <=
  deep or k-i <= deep.  Larger values of deep will usually yield
  shorter vectors, but the running increases exponentially.  

  NOTE: use of "deep" is obsolete, and has been "deprecated".
  It is recommended to use BKZ_FP to achieve higher-quality reductions.
  Moreover, the Givens versions do not support "deep", and setting
  deep != 0 will raise an error in this case.

* The optional parameter "check" is a function that is invoked after
  each size reduction with the current row as an argument.  If this
  function returns a non-zero value, the LLL procedure is immediately
  terminated.  Note that it is possible that some linear dependencies
  remain undiscovered, so that the calculated rank value is in fact
  too large.  In any case, zero rows discovered by the algorithm
  will be placed at the beginning, as usual.

  The check argument (if not zero) should be a routine taking
  a const vec_ZZ& as an argument and return value of type long.
  LLLCheckFct is defined via a typedef as:

     typedef long (*LLLCheckFct)(const vec_ZZ&);

  See the file subset.c for an example of the use of this feature.

* The optional parameter "verbose" can be set to see all kinds of fun
  things printed while the routine is executing.  A status report is
  printed every once in a while, and the current basis is optionally
  dumped to a file.  The behavior can be controlled with these global
  variables:

     extern char *LLLDumpFile;  // file to dump basis, 0 => no dump; 
                                // initially 0

     extern double LLLStatusInterval; // seconds between status reports 
                                      // initially 900s = 15min



 
************* Calling Syntax for BKZ routines ***************

long 
[G_]BKZ_{FP,QP,QP1,XD,RR} (mat_ZZ& B, [ mat_ZZ& U, ] double delta=0.99,
                          long BlockSize=10, long prune=0, 
                          LLLCheckFct check = 0, long verbose = 0);

These functions are equivalent to the LLL routines above,
except that Block Korkin-Zolotarev reduction is applied.
We describe here only the differences in the calling syntax.

* The optional parameter "BlockSize" specifies the size of the blocks
  in the reduction.  High values yield shorter vectors, but the
  running time increases exponentially with BlockSize.
  BlockSize should be between 2 and the number of rows of B.

* The optional parameter "prune" can be set to any positive number to
  invoke the Volume Heuristic from [Schnorr and Horner, Eurocrypt
  '95].  This can significantly reduce the running time, and hence
  allow much bigger block size, but the quality of the reduction is
  of course not as good in general.  Higher values of prune mean
  better quality, and slower running time.  
  When prune == 0, pruning is disabled.
  Recommended usage: for BlockSize >= 30, set 10 <= prune <= 15.

* The QP1 variant uses quad_float precision to compute Gramm-Schmidt,
  but uses double precision in the search phase of the block reduction
  algorithm.  This seems adequate for most purposes, and is faster
  than QP, which uses quad_float precision uniformly throughout.


******************** How to choose? *********************

I think it is safe to say that nobody really understands
how the LLL algorithm works.  The theoretical analyses are a long way
from describing what "really" happens in practice.  Choosing the best
variant for a certain application ultimately is a matter of trial
and error.

The first thing to try is LLL_FP.
It is the fastest of the routines, and is adequate for many applications.

If there are precision problems, you will most likely get
a warning message, something like "warning--relaxing reduction".
If there are overflow problems, you should get an error message
saying that the numbers are too big.

If either of these happens, the next thing to try is G_LLL_FP,
which uses the somewhat slower, but more stable, Givens rotations.
This approach also has the nice property that the numbers remain
smaller, so there is less chance of an overflow.

If you are still having precision problems with G_LLL_FP,
try LLL_QP or G_LLL_QP, which uses quadratic precision.

If you are still having overflow problems, try LLL_XD or G_LLL_XD.

I haven't yet come across a case where one *really* needs the
extra precision available in the RR variants.

All of the above discussion applies to the BKZ variants as well.
In addition, if you have a matrix with really big entries, you might try 
using G_LLL_FP or LLL_XD first to reduce the sizes of the numbers,
before running one of the BKZ variants.

Also, one shouldn't rule out using the "all integer" LLL routines.
For some highly structured matrices, this is not necessarily
much worse than some of the floating point versions, and can
under certain circumstances even be better.


******************** Implementation notes *********************

For all the floating point variants, I use a "relaxed" size reduction
condition.  Normally in LLL one makes all |\mu_{i,j}| <= 1/2.
However, this can easily lead to infinite loops in floating point arithemetic.
So I use the condition |\mu_{i,j}| <= 1/2 + fudge, where fudge is 
a very small number.  Even with this, one can fall into an infinite loop.
To handle this situation, I added some logic that detects, at quite low cost,
when an infinite loop has been entered.  When that happens, fudge
is replaced by fudge*2, and a warning message "relaxing reduction condition"
is printed.   We may do this relaxation several times.
If fudge gets too big, we give up and abort, except that 
LLL_FP and BKZ_FP make one last attempt to recover:  they try to compute the
Gramm-Schmidt coefficients using RR and continue.  As described above,
if you run into these problems, which you'll see in the error/warning
messages, it is more effective to use the QP and/or Givens variants.

For the Gramm-Schmidt orthogonalization, lots of "bookeeping" is done
to avoid computing the same thing twice.

For the Givens orthogonalization, we cannot do so many bookeeping tricks.
Instead, we "cache" a certain amount of information, which
allows us to avoid computing certain things over and over again.

There are many other hacks and tricks to speed things up even further.
For example, if the matrix elements are small enough to fit in
double precision floating point, the algorithms avoid almost
all big integer arithmetic.  This is done in a dynamic, on-line
fashion, so even if the numbers start out big, whenever they
get small, we automatically switch to floating point arithmetic.

\**************************************************************************/




/**************************************************************************\

                         Other Stuff

\**************************************************************************/



void ComputeGS(const mat_ZZ& B, mat_RR& mu, vec_RR& c);

// Computes Gramm-Schmidt data for B.  Assumes B is an m x n matrix of
// rank m.  Let if { B^*(i) } is the othogonal basis, then c(i) =
// |B^*(i)|^2, and B^*(i) = B(i) - \sum_{j=1}^{i-1} mu(i,j) B^*(j).

void NearVector(vec_ZZ& w, const mat_ZZ& B, const vec_ZZ& a);

// Computes a vector w that is an approximation to the closest vector
// in the lattice spanned by B to a, using the "closest plane"
// algorithm from [Babai, Combinatorica 6:1-13, 1986].  B must be a
// square matrix, and it is assumed that B is already LLL or BKZ
// reduced (the better the reduction the better the approximation).
// Note that arithmetic in RR is used with the current value of
// RR::precision().

// NOTE: Both of these routines use classical Gramm-Schmidt
// orthogonalization.