(*:Name: NumberTheory`ShortestVector` *)

(*:Summary:  ShortestVector uses my variant of the Fincke-Pohst      *)
(*           algorithm to find the shortest vector in a lattice.     *)
(*           It requires NumberTheory`LLLalgorithm`, preferably      *)
(*           Version 1.8 or higher.                                  *)

(*:Author:        Wilberd van der Kallen, June 1996, November 1999   *)
(*                http://www.math.uu.nl/people/vdkallen/kallen.html  *)
(*                bugs to:   vdkallen@math.uu.nl                     *)

(*:Mathematica Version: 2.2 or 3.0 *)
(*:Package Version: 1.8.1 *)

(*:Sources:                                                          *)
(* Henri Cohen, A course in computational Algebraic Number Theory,   *)
(* Graduate Texts in Mathematics 138, Springer 1993.                 *)

(*:Discussion:                                                       *)
(* Given sufficient Accuracy, the algorithm finds a shortest vector  *)
(* although not necessarily in polynomial time. As a rule of thumb   *)
(* the WorkingPrecision should be a little more than is sufficient   *)
(* to distinguish the square length of the shortest vector from      *)
(* other integers. If one chooses WorkingPrecision larger than       *)
(* $MachinePrecision, or if one chooses WorkingPrecision->Automatic, *)
(* then the algorithm checks if the Accuracy stays in a safe range   *)
(* and a warning is given if it does not.                            *)
(* If one has a guess of a shortest vector, it may help to prepend   *)
(* that guess to the generating set. One may also wish to preprocess *)
(* with LatticeReduce.                                               *)

BeginPackage["ShortestVector`","LLLalgorithm`"]

ShortestVector::usage="If {b1,...,bn} is a generating set of
an integral lattice, then
ShortestVector[{b1,...,bn}] is
a shortest non-zero vector in the lattice."

(*:Examples:

In[1]:= <<ShortestVector.m

In[2]:= generatingset={{1,2,3},{4,5,6},{7,8,9},{10,11,12}};

In[3]:= ShortestVector[generatingset]

Out[3]= {-2, -1, 0}

In[4]:= ShortestVector[10^25 generatingset]

Out[4]= {-20000000000000000000000000, -10000000000000000000000000, 0}
  (* No warning is given when working with machine precision. *)

In[5]:= Check[ShortestVector[10^25 generatingset,
               WorkingPrecision->3 $MachinePrecision],$Failed]

ShortestVector::prec:
   WorkingPrecision of 48 may not suffice to find true shortest vector.

Out[5]= $Failed

In[6]:= Check[ShortestVector[10^25 generatingset,
               WorkingPrecision->4 $MachinePrecision],$Failed]

Out[6]= {-20000000000000000000000000, -10000000000000000000000000, 0}

In[7]:= Check[ShortestVector[10^25 generatingset,
               WorkingPrecision->Automatic],$Failed]

Out[7]= {-20000000000000000000000000, -10000000000000000000000000, 0}

*)

Unprotect[ShortestVector];

Options[ShortestVector]={
RatioLLLCondition->{3/4,19/20,99/100},
WorkingPrecision->$MachinePrecision};


Begin["ShortestVector`Private`"]

(* Algorithm 2.7.5 ("Short Vectors" after Fincke & Pohst) from H.Cohen's book *)

shortVectors::usage:="Let C be a constant and Q the positive definite
quadratic form
Q[x_]=Sum[q[[i,i]](x[[i]]+Sum[q[[i,j]]x[[j]],{j,i+1,n}])^2,{i,n}].
Then shortVectors[q,C] is a list of nonzero vectors with integer coefficients
and Q-value smaller than about C, listing only one from each pair x, -x.
This algorithm is only one step in the approach of Fincke & Pohst and
should be tried only for very small cases or after preprocessing with a
reduction procedure."

shrinking::usage:="shrinking is an option in shortVectors and
shortVectorsGram. shrinking->True makes that a vector is
skipped when it is longer than the preceding one."

Options[shortVectors]={WorkingPrecision->$MachinePrecision,shrinking->False};

shortVectors[quadraticdata_,bound_,options___Rule]:=
Module[ 
        {q,c,n,i,j,t,u,l,x,z,sol={},workingPrec,highprec,shrink,
        computebounds,mainloop,solutionfound,next,finished
        }
        ,
        shrink=True===(shrinking/.{options}/.Options[shortVectors]);
        workingPrec=(WorkingPrecision/.{options}/.Options[shortVectors]);
        highprec=(workingPrec>$MachinePrecision);
	q=N[quadraticdata,workingPrec];
	c=N[Min[1001/1000bound,bound+1/2],workingPrec];
	If[c<=bound,c=N[1001/1000bound,workingPrec]];
        i=n=Length[q];
        u=l=x=t=Range[n];
        t[[i]]=c;
        u[[i]]=0;
	computebounds:=(
                       z=Sqrt[t[[i]]/q[[i,i]]];
		       If[highprec&&(2Accuracy[z+u[[i]]]<
			    4+  (* 4 is on the safe side *)
                            Last[MantissaExponent[q[[i,i]]t[[i]]]])
                         ,Message[ShortestVector::prec,workingPrec];
			  highprec=False
                         ];(* This assumes integer valued form *)
                       l[[i]]=Floor[z-u[[i]]];
                       x[[i]]=Ceiling[-z-u[[i]]]-1;
                       next:=mainloop;
                       );
        mainloop:=(
                  x[[i]]=x[[i]]+1;
                  If[x[[i]]>l[[i]],
                    i++;next:=mainloop,
                    If[i>1,
                      t[[i-1]]=t[[i]]-q[[i,i]](x[[i]]+u[[i]])^2;
                      i--;
                      u[[i]]=Sum[q[[i,j]]x[[j]],{j,i+1,n}];
                      next:=computebounds;
                      ,
                      next:=solutionfound
                      ];
                    ];
                  );
	solutionfound:=(
			If[x==0x,
			   finished=True,
			   If[shrink,
			      t=t-(t[[1]]-q[[1,1]](x[[1]]+u[[1]])^2)
			     ];
			   sol=Append[sol,x];
			   next:=mainloop
			  ]
		       );
	For[finished=False;next:=computebounds;,!finished,Null,next];
	computebounds=mainloop=solutionfound=Null;(* free Temporary Symbols *)
        sol
        ];



ShortestVector[m___]:=ShortestVectorAnswer[m]/;ShortestVectorTest[m]

ShortestVectorAnswer[m_,options___Rule]:=
Module[{basis,workingPrec,opts},
	opts=Join[{options},Options[ShortestVector]];
	basis=First[ExtendedLatticeReduce[m,opts]];
        workingPrec=(WorkingPrecision/.opts);
	If[workingPrec===Automatic,
	   workingPrec=$MachinePrecision+
                       Last[MantissaExponent[1.*Max[(#.#&)/@(basis)]]];
	   opts=Join[{WorkingPrecision->workingPrec},opts];
	  ];
	Last[shortVectorsGram[
			N[basis.Transpose[basis],workingPrec],
			N[Min[(#.#&)/@(basis)],workingPrec],
			shrinking->True,opts
			     ].basis
            ]
      ];

ShortestVectorAnswer[m__,{rules___},more___]:=
ShortestVectorAnswer[m,rules,more]

ShortestVectorTest[m_,options___Rule]:=
Module[{basis,failed,error,workingPrec},
	workingPrec=(WorkingPrecision/.{options}/.Options[ShortestVector]);
	failed=False;
	(* begin check argument *)
	If[!((IntegerQ[workingPrec])||(workingPrec===Automatic)),
		Message[ShortestVector::val,workingPrec];failed=True
	  ];
	basis=m;
	error:=(Message[ShortestVector::argint,basis];failed=True;error=Null);
	If[MatrixQ[basis]
	   ,
           If[(Length[basis[[1]]])===0,error];
   	   Map[(If[ #[[0]]===Integer,Null,error])&,basis,{2}]
  	   ,
  	   error
  	  ];
	(* end check argument *)
        !failed
      ];

ShortestVectorTest[m__,{rules___},more___]:=
ShortestVectorTest[m,rules,more]

ShortestVectorTest[___]:=(Message[ShortestVector::argw];False)

ShortestVector::argint:=
"Argument `1` of ShortestVector should be a matrix with integral
entries."

ShortestVector::argw:="ShortestVector called with wrong
number of arguments. One argument is expected, not counting options.";

ShortestVector::val:="`1` is not a valid value for the option
WorkingPrecision.";

ShortestVector::prec:=
"WorkingPrecision of `1` may not suffice to find true shortest vector."

shortVectorsGram::usage="If {b1,...,bn} is a generating set of
an integral lattice, and gram={b1,...,bn}.Transpose[{b1,...,bn}]
is its Gram matrix, then
shortVectorsGram[gram,C].{b1,...,bn} is
a list of half the non-zero vectors in the lattice shorter than roughly
Sqrt[C].
This algorithm is only one step in the approach of Fincke & Pohst and
should be tried only for very small cases or after preprocessing with a
reduction procedure."

shortVectorsGram[m_,bound_,options___Rule]:=(* m=basis.Transpose[basis] *)
Module[
	{b,bb,i,j,n,gram=m,mu,bst}
	,
	n=Length[gram];
	bst=b=IdentityMatrix[n];
	bb=Range[n];
	mu=Table[0,{i,n},{j,n}];
	Do[
	  Do[
	     mu[[i,j]]=(b[[i]].gram.bst[[j]])/bb[[j]];
	     bst[[i]]=bst[[i]]-mu[[i,j]] bst[[j]];
	     ,{j,i-1}
	    ];
	  mu[[i,i]]=bb[[i]]=bst[[i]].gram.bst[[i]];
	  ,{i,n}
	  ];
	shortVectors[Transpose[mu],bound,options]
      ];

shortVectorsGram[m__,{rules___},more___]:=
shortVectorsGram[m,rules,more]

End[]

Protect[ShortestVector];

EndPackage[]


Null
