(*: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]:= <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