NTL库之LLL，BKZ约化算法

可选函数类型

本文件中包含格基约化算法所需要的必备函数，这些功能包含定点数版本（exact arithmetic variants，特点是虽然慢不过是精确解）以及浮点数版本（floating point variants，快，但是近似解）如果读者想要进一步了解LLL算法，请移步[H. Cohen, A Course in Computational Algebraic Number Theory, Springer, 1993].LLL算法最早是在文献 [A. K. Lenstra, H. W. Lenstra, andL. Lovasz, Math. Ann. 261 (1982), 515-534].中提出的。

可选函数类型

long LLL(ZZ& det2, mat_ZZ& B, long verbose = 0);
long LLL(ZZ& det2, mat_ZZ& B, mat_ZZ& U, long verbose = 0);
long LLL(ZZ& det2, mat_ZZ& B, long a, long b, long verbose = 0);
long LLL(ZZ& det2, mat_ZZ& B, mat_ZZ& U, long a, long b, long verbose = 0);

B是一个m*n的矩阵，我们将其看成是m个行向量组成，每个行向量分量个数都是n。 行数m小于等于n。并且行向量之间并不需要线性无关。B被转化成一个LLL约化基，并且返回的值是B的秩，我们将其记做r。

更具体的来说，我们是在B上进行行变换，新B的非零行是原始的B矩阵的LLL约化基。默认的规约参数为delta = 3/4,这个值意味着第一个非零的基向量的长度不超过格中最短向量的2**r-1次。（r是矩阵的秩）Det2指的是格的判别式（determinant）的平方，当r=n（列满秩情况下）sqrt(det2)通常是一个整数。

在第二个版本中，U被设定为转换矩阵，满足U是一个单模的mm矩阵，Uold-B=new B,我们得注意到U的前m-r行是old-B的kernel的基。

第三个和第四个版本允许对任意的delta=a/b，进行格基约化，其中1/4<a/b<=1，a和b 均为正整数。使用参数delta约化得到的一组基，其第一个非零向量的平方长度不会超过格上最短向量的1/(delta-1/4)^(r-1)倍。 （我们下面给出的Schnorr和Euchner的文章中会体现这一点）
这里采用的算法本质上还是Cohen的书里面给出的算法。

// Some variations:

long LLL_plus(vec_ZZ& D, mat_ZZ& B, long verbose = 0);
long LLL_plus(vec_ZZ& D, mat_ZZ& B, mat_ZZ& U, long verbose = 0);
long LLL_plus(vec_ZZ& D, mat_ZZ& B, long a, long b, long verbose = 0);
long LLL_plus(vec_ZZ& D, mat_ZZ& B, mat_ZZ& U, long a, long b,long verbose = 0);

上面给出的规约函数版本会多返回1bit关于约化基的信息。如果r是B的秩，则D就是一个长度为r+1的向量，其中D[0]=1,并且对于i=1,…,r，D[i]/D[i-1]等于第i个对应于LLL约化基B的非零行的Gram-Schmidt变换之后长度的平方（拿这里做一下实验）。具体来说，D[r]等于LLL格基约化之后得到的det2。

计算图像和kernel

long image(ZZ& det2, mat_ZZ& B, long verbose = 0);
long image(ZZ& det2, mat_ZZ& B, mat_ZZ& U, long verbose = 0);

该函数计算了B的图像，使用的是LLL的简化版本，在本函数中，使用了常规的“尺寸缩小”技术，不过仅当发现了线性相关的对的时候才会交换向量，我暂时没有在文献中发现这种处理方法，但是这种方法在实际应用中表现非常好，并且可以在合理时间内控制循环次数。和上面的LLL实现一样，被返回的值是B的秩r，并且前面的m-r行是0.U是一个单模的mm矩阵，满足Uold-B= new-B. det2的含义和上面一致。

注意到U的前m-r中构成了old-B的kernel的一组基。这是一个非常实用的计算kernel的算法。人们同样可以给kernel使用image()函数以便获得在某种意义上更短的kernel的基（线性无关，但无论如何能把基的尺寸缩小一点点）。即便是对于更短的kernel的基向量，人们也可以使用LLL()。

在格中寻找一个向量

long LatticeSolve(vec_ZZ& x, const mat_ZZ& A, const vec_ZZ& y, long reduce=0);

LatticeSolve函数，传入参数为整形向量x，整形矩阵A，整形向量y，以及长整型reduce。

该函数测试了如果给定A和y，是否存在x满足xA=y。 如果满足，则x就是一个解，返回1，否则x还是传入的值，返回值为0.

最优参数降低了对解向量的“质量”的控制，如果A的行是线性相关的，则存在很多组解。降低控制被转移到了寻找“较短的”x的问题。

Reduce = 0: 为寻找一个短的解不需要特定的处理

Reduce = 1： 需要在kernel(A)上运行一个简单的“size reduction”，这种处理的效率很高，并且默认可能得到一个比较短的解，但是不一定非常接近最优解。

Reduce = 2：在kernel(A)上运行LLL算法，这可能比起其他选项来说运行的慢一点，但是得到的解被证明非常接近最优解。 更加精确的说，如果kernel(A)的秩为k，那么得到的解的平方不会超过最优解的max(1,2**(k-2))倍。 该算法内部使用了Babai的近似最近向量算法。

当然，如果A的行是线性无关的，则reduce使用什么值就无所谓了，因为在这个时候，得到的解都是唯一解。

我们得注意到无论reduce是什么值，算法运行时间都是多项式级别的，并且解向量的bit长度是输入长度的多项式级别。

浮点数变体

存在LLL算法的变体：可以选择精度，正交手段以及约化条件。

大量的选择可能看起来有点让人迷惑，所以请阅读下面的参数选择指导：

精度：

FP：双精度浮点(double)类型
QP：四倍精度浮点（类,quad_float,quasi quadruple） 当四舍五入导致的误差可能会导致严重问题的时候可用
XD:xdouble，扩展指数双精度浮点（extended exponent doubles），当数据过大的时候可以使用
RR：RR（任意精度的浮点数类型） 对于较大的精度和量级会比较好用

一般来说，FP的选择会是最快的，但是可能会导致四舍五入的误差以及数据溢出。

正交化策略：

-经典的Gram-Schmidt正交化，这种选择使用经典的计算Gram-Schimidt正交化计算方法，这种计算方法非常快，但是可能导致计算稳定性问题。这种方法最早是Schnorr和Euchner提出的。详见：[C. P. Schnorr and M. Euchner, Proc. Fundamentals of Computation Theory, LNCS 529, pp. 68-85, 1991]. 这里实现的版本和最早的版本有本质区别，改版之后提升了数据的稳定性以及实现性能。

-Givens正交化：这种正交化运行稍微慢一点，但是一般来说运行的更稳定一些，并且人们更加偏爱使用这种算法，如果想要更加详细了解这种算法， 请参考文献 [G. Golub and C. van Loan, Matrix Computations, 3rd edition,Johns Hopkins Univ. Press, 1996]的第五章。

约化条件：

—LLL：经典的LLL约化条件

—BKZ：Block Korkin Zolotarev规约，这种算法更慢一点，不过可以得到更加高质量的基，也就是说，能够得到更短的向量。 可以看一下Schonorr-Euchner的文章中关于本算法的描述，该算法本质上扩展了LLL约减条件，将分块数为2拓展到了更大的分块数。

LLL调用语法

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);

算法名分为LLL_FP,LLL-QP,LLL_XD,LLL_RR,分别对应不同的计算精度，传入的参数分别为整形矩阵B，可选项为整形矩阵U，根hermite因子（root hermite factor）delta，长整型deep，LLLCheckFct类型参数check，长整型参数verbose。

*[…]表达式是可选项，{…}指的是在多个可选项中选择其中一个
*返回值是B的rank（如果check不是0则继续往下看）
*可选前缀G_表示使用Givens旋转，否则，使用经典的Gram-Schmidt正交化算法。
*FP,QP,XD,RR决定了使用的数据精度
*如果给定了可选参数U，则U被用作转换矩阵

 U * old_B = new_B

*可选参数delta是约化参数，其设定区间是[0.5,1)，这个值越接近1，则得到的向量越短，并且同样提升了计算稳定性，不过运行时间增加了，我们一般推荐取delta = 0.99

*可选参数“deep”可以被设定为任意正整数，该正整数允许人们“更深”的将第k行插入第i行，前提条件是i<=deep或者k-i<=deep。更大的deep通常会得出更短的向量，不过运行时间会指数级别上升。
注意：使用“deep”实际有点过时了，并且已经被弃用了，我们还是推荐使用BKZ_FP以便达到更高质量的格基约化。此外，Givens版本并不支持deep，如果设置deep为非零的值会抛出异常。
*可选项check在每次尺寸约减之后调用，参数为当前的行。 如果这个函数返回一个非零的值，LLL过程直接终止。 我们得注意到有可能还有一些线性相关关系没有被发现，从而计算出来的rank可能过大。 在任意情况下，算法发现的零行会被放在最开始。

Check参数（如果非零）应该是一个过程：该过程的输入是一个常指针vec_ZZ&，作为一个参数，并且返回值是长整形。LLLCheckFct通过typedef定义为：

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

详情请见subset.cpp，给出了使用这个特征的例子。

*设定可选参数“verbose”（冗长？），从而能够看到各种运行过程中输出的有意思的东西。 while循环每执行一次都会输出状态报告，并且会将当前的基输出到文件中，可以通过下面的全局变量

 extern thread_local char *LLLDumpFile;  
 //说明输出格基的文件位置，如果指定为0则认为不向中间文件输出格基
 extern thread_local double LLLStatusInterval; 
 // 产生两个中间报告的时间间隔，初始值为900s=15min

BKZ 调用语法

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);

5类调用接口，分别为BKZ_FP，BKZ_QP,BKZ_QP1,BKZ_XD,BKZ_RR，输入参数分别为整形矩阵指针B，整形矩阵指针U（转换矩阵），root hermite factor 0.99，分块数BlockSize，默认值为10，长整形 prune，默认为0，LLLCheckFct类参数check=0，长整型verbose=0.

这些函数内部使用BKZ约化算法。我们在下面仅仅描述调用语法的区别。

*可选参数“BlockSize”说明了规约中使用的分块大小，分块数越大，得到的向量越短，但是运行时间会随着分块数增长而指数级别增加。BlockSize的取值应该在2和B的行数之间选取。

*可选参数prune（剪枝）：可以被设定为任意正数，以便调用体积启发式算法，详情请见schnorr和Horn
er在EuroCrypt95年的文章。设置这个值能够在很大程度上降低运行时间，从而能够允许人们选择更大的值，不过格基约化的质量比起不使用剪枝来说要差一点。 剪枝参数prune的值越高，意味着更高的质量，运行时间更长，如果prune=0，禁用剪枝。
推荐剪枝参数（prune）的值： 当BlockSize>=30,设定10<=prune<=15

*QP1变体会使用quad_float精度（四倍精度浮点数）以便计算Gram-Schmidt参数，但是在块内约化的搜索阶段使用了双精度浮点数。这种精度对于大部分情况都足够了，并且比起QP来说要快，因为QP中整体都使用四倍精度。

如何选择？

我认为这么说比较稳妥：没有人真的理解LLL算法怎么工作的。理论分析距离真正去描述运行还远的很。 为特定的应用场景选择合适的变体归根结底还是一个试错的过程。

第一个值得尝试的事情是LLL_FP，这是所有可选项中最快的部分， 并且对很多场景来说足够了。

如果存在精度问题，我们可能会得到一个警告信息，警告信息可能是“warning–relaxing reduction”。如果存在溢出问题，则可能得到错误信息，错误报告为数据过大。

如果上方两种情况都没有发生，下面需要尝试的事情是G_LLL_FP,该变体使用的Givens
旋转。这种方法同样也有非常好的性质：计算过程中的数都很小，从而出现数据移除的情况不多。
如果在G_LLL_FP中，仍然存在数据精度问题，可以尝试使用LLL_QP,G_LLL_QP,这些变体中使用了四倍精度。

如果您仍然存在数据溢出问题，请尝试使用LLL_XD或者G_LLL_XD.
我目前还没有遇到人们“真的”需要在RR版本中使用更高的精度。

上面的全部讨论对于BKZ变体来说都同样适用，此外，如果你手中拥有一个分量真的非常大的矩阵，那你可能需要使用G_LLL_FP,或者LLL_XD,以便首先降低数据的尺寸，随后运行BKZ算法中的一个。

同样，人们不能排除使用全部为整数的“LLL”过程，对于某些高度结构化的矩阵来说，这些变体不一定比浮点数版本差，并且在某些特定情况下甚至更好一点。

实现中的注意事项

对于全部的浮点数变体，我使用了一个更加“放宽”的尺寸降低条件。正常来说在LLL中，人们会令.但是，这个条件很容易在浮点数版本中导致死循环，所以所以我将这个条件进行了放宽，使用，其中参数fudge是一个非常小的数。

即便使用了放宽的条件，也还是有可能会陷入无限循环中，为了解决这个情况，我们会将fudge放大到fudge*2，并且返回一个警告信息“relaxing reduction condition”。

我们可能会多次进行上面的条件放宽过程。如果fudge变得太大，我们会放弃这个过程并且中止，除了LLL_FP和BKZ_FP会进行最后一次尝试以便恢复：这两个变体尝试使用RR数据类型计算Gram-Schmidt系数，和上面描述的一样，如果你遇到了这些问题（在error中或者Warning中能看到）中的一种，使用QP变体或者Givens变体可能更加高效。

对于Gram-Schmidt正交化，我们会进行很多非常有意思的中间记录，以便避免重复计算。

对于Givens正交化，我们不能做这么多中间记录，所以，我们会“缓存”特定两的信息，这些信息允许我们避免重复计算。

存在很多其他的解决方案和其他技巧以便进一步加速计算过程。举例来说，如果矩阵中的蒜素足够小，能够放进双倍精度浮点数，那么算法基本能够避免全部的大整数计算，这种过程是动态的、在线的，从而即便是当起始数据非常大的时候，只要在数据处理过程中这些数据开始变小，我们就会自动切换到浮点数据。

其他

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

// ComputeGS，参数为ZZ类型的B，以及RR类型的投影参数mu，RR类型的向量，记录了约化之后向量长度的平方。
//ComputeGS会为B计算Gram-Schmidt数据，假定B是一个mn的矩阵，秩为m，我们令{b^{*(i)}是正交基，那么得到的c(i)=|B}(i)|^2,并且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);

NearVector的参数为ZZ型的向量w，以及ZZ型的矩阵B，ZZ型的向量a，该函数计算最接近B张成的格a的最近向量，使用Babai等人提出的最近平面算法。B必须是一个平方矩阵，并且假定B已经是LLL或者BKZ约化基（约化程度越高，近似程度越好）我们首先注意到RR中的计算的精度为RR:precision()中的精度。

注意：上面两个过程中使用的都是Gram-Schmidt正交化。