LLL算法的C++实现

LLL.h

#pragma once

#include "tools.h"

#include <iostream>
#include <vector>

using namespace std;


/*
	基向量
*/
typedef vector<double> Vec;

/*
	打印行向量
*/
ostream& operator<<(ostream& cout, Vec& v);

/*
	格基
*/
typedef vector<Vec> Base;

/*
	按照行向量格式打印格基
*/
ostream& operator<<(ostream& cout, Base& B);

/*
	格基初始化
*/
void Init(Base& B, int32 n);


/*
	dst = v1 + v2
*/
void Add(Vec& dst, Vec& v1, Vec& v2, int32 n);


/*
	dst = v1 - v2
*/
void Sub(Vec& dst, Vec& v1, Vec& v2, int32 n);


/*
	dst = scale * v
*/
void Mul(Vec& dst, double scale, Vec& v, int32 n);


/*
	dst = v1 + scale * v2
*/
void Add(Vec& dst, Vec& v1, double scale, Vec& v2, int32 n);


/*
	dst = v1 - scale * v2
*/
void Sub(Vec& dst, Vec& v1, double scale, Vec& v2, int32 n);


/*
	向量内积
*/
double Inner(Vec& v1, Vec& v2, int32 n);


/*
	v1在v2上投影
*/
double Proj(Vec& v1, Vec& v2, int32 n);


/*
	Gram-Schmidt正交化
*/
void GramSchmidt(Base& dst, Base& src, int32 n);


/*
	Hadamard比率
*/
double Hadamard(Base& B, int32 n);


/*
	LLL算法
*/
void LLL(Base& dst, Base& src, int32 n, double delta = 0.75);


/*
	判断是否正交
*/
bool IsOrthogonal(Base& B, int32 n);

/*
	格基正交性
*/
void BaseOrthogonal(Base& B, int32 n);


/*
	标准化
*/
void Stand(Base& dst, Base& src, int32 n);


/*
	转化为劣质基
*/
void ToBadBase(Base& dst, Base& src, int32 n, int32 bnd = 1000);

LLL.cpp

#include "LLL.h"


/*
	误差
*/
#define eps 1e-10


/*
	约等于
*/
#define approx(a,b) (abs((a)-(b))<eps)


/*
	打印行向量
*/
ostream& operator<<(ostream& cout, Vec& v)
{
	int32 n = v.size();
	double* p = v.data();

	printf("[");
	for (int32 i = 0; i < n; i++)
	{
		printf(" %+15.3lf ", *p);
		p++;
	}
	printf("]\t|\t%+15.3lf", sqrt(Inner(v, v, n)));

	return cout;
}


/*
	按照行向量格式打印格基
*/
ostream& operator<<(ostream& cout, Base& B)
{
	int32 n = B.size();

	printf("[");
	for (int32 i = 0; i < n; i++)
	{
		printf("\n\t");
		cout << B[i];
	}
	printf("\n]\n");

	return cout;
}



/*
	格基初始化
*/
void Init(Base& B, int32 n)
{
	B.resize(n);
	for (int32 i = 0; i < n; i++)
	{
		B[i].resize(n);
		memset(B[i].data(), 0, sizeof(double)*n);
	}
}


/*
	dst = v1 + v2
*/
void Add(Vec& dst, Vec& v1, Vec& v2, int32 n)
{
	dst.resize(n);
	double* p1 = v1.data();
	double* p2 = v2.data();
	double* p3 = dst.data();

	for (int32 i = 0; i < n; i++)
	{
		*p3 = *p1 + *p2;
		p1++;
		p2++;
		p3++;
	}
}


/*
	dst = v1 - v2
*/
void Sub(Vec& dst, Vec& v1, Vec& v2, int32 n)
{
	dst.resize(n);
	double* p1 = v1.data();
	double* p2 = v2.data();
	double* p3 = dst.data();

	for (int32 i = 0; i < n; i++)
	{
		*p3 = *p1 - *p2;
		p1++;
		p2++;
		p3++;
	}
}


/*
	dst = scale * v
*/
void Mul(Vec& dst, double scale, Vec& v, int32 n)
{
	dst.resize(n);
	double* p1 = v.data();
	double* p2 = dst.data();

	for (int32 i = 0; i < n; i++)
	{
		*p2 = scale * *p1;
		p1++;
		p2++;
	}
}


/*
	dst = v1 + scale * v2
*/
void Add(Vec& dst, Vec& v1, double scale, Vec& v2, int32 n)
{
	dst.resize(n);
	double* p1 = v1.data();
	double* p2 = v2.data();
	double* p3 = dst.data();

	for (int32 i = 0; i < n; i++)
	{
		*p3 = *p1 + scale * *p2;
		p1++;
		p2++;
		p3++;
	}
}


/*
	dst = v1 - scale * v2
*/
void Sub(Vec& dst, Vec& v1, double scale, Vec& v2, int32 n)
{
	dst.resize(n);
	double* p1 = v1.data();
	double* p2 = v2.data();
	double* p3 = dst.data();

	for (int32 i = 0; i < n; i++)
	{
		*p3 = *p1 - scale * *p2;
		p1++;
		p2++;
		p3++;
	}
}


/*
	向量内积
*/
double Inner(Vec& v1, Vec& v2, int32 n)
{
	double* p1 = v1.data();
	double* p2 = v2.data();

	double res = 0;
	for (int32 i = 0; i < n; i++)
	{
		res += *p1 * *p2;
		p1++;
		p2++;
	}

	return res;
}


/*
	v1在v2上投影
*/
double Proj(Vec& v1, Vec& v2, int32 n)
{
	return Inner(v1, v2, n) / Inner(v2, v2, n);
}


/*
	Gram-Schmidt正交化
*/
void GramSchmidt(Base& dst, Base& src, int32 n)
{
	if (dst != src)
	{
		Init(dst, n);
	}

	for (int32 i = 0; i < n; i++)
	{
		dst[i] = src[i];
		for (int32 j = 0; j < i; j++)
		{
			double mu = Proj(dst[i], dst[j], n);
			Sub(dst[i], dst[i], mu, dst[j], n);
		}
	}
}


/*
	Hadamard比率
*/
double Hadamard(Base& B, int32 n)
{
	Base BB;
	GramSchmidt(BB, B, n);
	double r = 1, s = 1;
	for (int32 i = 0; i < n; i++)
		r *= sqrt(Inner(BB[i], BB[i], n));
	for (int32 i = 0; i < n; i++)
		s *= sqrt(Inner(B[i], B[i], n));
	return r / s;
}


/*
	LLL算法
*/
void LLL(Base& dst, Base& src, int32 n, double delta)
{
	if (dst != src)
	{
		Init(dst, n);
		for (int32 i = 0; i < n; i++)
			dst[i] = src[i];
	}

	Base B;
	Init(B, n);

	int32 k = 1;
	Vec tmp(n);

	while (k < n)
	{
		if (k == 1)
			B[0] = dst[0];

		// Size条件
		for (int32 j = k - 1; j >= 0; j--)
		{
			double mu = round(Proj(dst[k], B[j], n));
			Sub(dst[k], dst[k], mu, dst[j], n);
		}

		// Gram - Schmidt正交化
		B[k] = dst[k];
		for (int32 j = 0; j < k; j++)
		{
			double mu = Proj(B[k], B[j], n);
			Sub(B[k], B[k], mu, B[j], n);
		}

		// Lovász条件
		double mu = Proj(dst[k], B[k-1], n);
		if (Inner(B[k], B[k], n) >= (delta - mu * mu) * Inner(B[k - 1], B[k - 1], n))
			k += 1;
		else
		{
			// 交换位置，使得范数下降
			tmp = dst[k];
			dst[k] = dst[k - 1];
			dst[k - 1] = tmp;

			k = k >= 2 ? k - 1 : 1;
		}
	}
}


/*
	判断是否正交
*/
bool IsOrthogonal(Base& B, int32 n)
{
	for (int32 i = 1; i < n; i++)
		for (int32 j = 0; j < i; j++)
			if (!approx(Inner(B[i], B[j], n), 0))
				return 0;
	return 1;
}


/*
	格基正交性
*/
void BaseOrthogonal(Base& B, int32 n)
{
	for (int32 i = 0; i < n; i++)
	{
		printf("[");
		for (int32 j = 0; j < n; j++)
			printf(" %+.5lf ", Inner(B[i], B[j], n)/(sqrt(Inner(B[i], B[i], n)) * sqrt(Inner(B[j], B[j], n))));
		printf("]\n");
	}
}


/*
	标准化
*/
void Stand(Base& dst, Base& src, int32 n)
{
	if (dst != src)
	{
		Init(dst, n);
	}

	for (int32 i = 0; i < n; i++)
	{
		double L2 = sqrt(Inner(src[i], src[i], n));
		Mul(dst[i], 1 / L2, src[i], n);
	}
}


/*
	转化为劣质基
*/
void ToBadBase(Base& dst, Base& src, int32 n, int32 bnd)
{
	Base B;
	vector<uint32> rd(n);
	
	Init(B, n);
	for (int32 i = 0; i < n; i++) //右乘上三角幺模矩阵
	{
		Add(B[i], B[i], src[i], n);

		Uniform_int32(rd, n);
		for (int32 j = 0; j < i; j++)
		{
			Add(B[i], B[i], rd[j] % bnd, src[j], n);
		}
	}

	Init(dst, n);
	for (int32 i = 0; i < n; i++) //右乘下三角幺模矩阵
	{
		Add(dst[i], dst[i], B[i], n);

		Uniform_int32(rd, n);
		for (int32 j = i+1; j < n; j++)
		{
			Add(dst[i], dst[i], rd[j] % bnd, B[j], n);
		}
	}
}

Test.cpp

#include "LLL.h"

int main()
{
	int32 n = 8;
	Base B1, B2, B3, B4, B5;

	Init(B1, n);
	for (int32 i = 0; i < n; i++)
		for (int32 j = 0; j < n; j++)
			B1[i][j] = Uniform_int32_dist(PRG) % 10;

	//B1 = { {1,0,0,0},{0,1,0,0},{0,0,1,0},{0,0,0,1} };

	cout << "优质基 = \n" << B1 << endl;
	cout << "Hadamard比率" << Hadamard(B1, n); pn; pn;

	/*GramSchmidt(B3, B1, n);
	cout << "Gram-Schmidt正交基 = \n" << B3 << endl;
	cout << "IsOrthogonal = " << IsOrthogonal(B3, n) << endl << endl;

	Stand(B4, B3, n);
	cout << "标准正交基 = \n" << B4 << endl;*/

	ToBadBase(B2, B1, n, 1000);
	cout << "劣质基 = \n" << B2 << endl;
	cout << "Hadamard比率" << Hadamard(B2, n); pn; pn;

	/*GramSchmidt(B3, B2, n);
	cout << "Gram-Schmidt正交基 = \n" << B3 << endl;
	cout << "IsOrthogonal = " << IsOrthogonal(B3, n) << endl << endl;

	Stand(B4, B3, n);
	cout << "标准正交基 = \n" << B4 << endl;*/

	double delta = 0.3;
	LLL(B5, B2, n, delta);
	cout << "delta=0.3 格基约化 = \n" << B5 << endl;
	cout << "Hadamard比率" << Hadamard(B5, n); pn; pn;
	cout << "正交性：\n"; BaseOrthogonal(B5, n); pn;

	delta = 0.75;
	LLL(B5, B2, n, delta);
	cout << "delta=0.75 格基约化 = \n" << B5 << endl;
	cout << "Hadamard比率" << Hadamard(B5, n); pn; pn;
	cout << "正交性：\n"; BaseOrthogonal(B5, n); pn;

	delta = 0.95;
	LLL(B5, B2, n, delta);
	cout << "delta=0.95 格基约化 = \n" << B5 << endl;
	cout << "Hadamard比率" << Hadamard(B5, n); pn; pn;
	cout << "正交性：\n"; BaseOrthogonal(B5, n); pn;

	return 0;
}

结果

优质基 =
[
        [          +2.000           +6.000           +4.000           +7.000           +4.000           +7.000           +6.000           +9.000 ]      |               +16.941
        [          +4.000           +8.000           +7.000           +4.000           +3.000           +5.000           +1.000           +3.000 ]      |               +13.748
        [          +3.000           +5.000           +2.000           +6.000           +2.000           +9.000           +2.000           +8.000 ]      |               +15.067
        [          +9.000           +1.000           +8.000           +7.000           +1.000           +9.000           +1.000           +3.000 ]      |               +16.941
        [          +3.000           +7.000           +5.000           +1.000           +7.000           +2.000           +7.000           +6.000 ]      |               +14.900
        [          +2.000           +0.000           +4.000           +4.000           +5.000           +6.000           +7.000           +5.000 ]      |               +13.077
        [          +2.000           +4.000           +3.000           +1.000           +6.000           +8.000           +2.000           +2.000 ]      |               +11.747
        [          +8.000           +9.000           +6.000           +8.000           +6.000           +4.000           +8.000           +2.000 ]      |               +19.105
]

Hadamard比率0.000111934

劣质基 =
[
        [   +29826647.000    +37125995.000    +37136826.000    +42874336.000    +26322723.000    +54094174.000    +30077435.000    +49320020.000 ]      |        +111541484.591
        [   +36591579.000    +46749101.000    +45688882.000    +49967509.000    +32755485.000    +65431678.000    +34826265.000    +57908153.000 ]      |        +134272473.699
        [   +27465374.000    +32943703.000    +34383356.000    +36178300.000    +25282952.000    +48071604.000    +26999086.000    +41918182.000 ]      |         +98827679.812
        [   +13624911.000    +14238206.000    +15766783.000    +17251544.000    +10607007.000    +22951626.000    +11442741.000    +19162811.000 ]      |         +45519306.969
        [    +8693029.000    +11060379.000    +10815147.000    +11879245.000     +8062262.000    +15368186.000     +8907054.000    +14225678.000 ]      |         +32231527.085
        [   +10145545.000    +10860657.000    +12611692.000    +13929246.000     +8897025.000    +17924641.000    +10266913.000    +15466200.000 ]      |         +36318415.121
        [    +1280627.000     +1843441.000     +1817097.000     +1540282.000     +1652868.000     +2207500.000     +1654010.000     +2071904.000 ]      |          +5034496.452
        [      +11738.000       +16953.000       +16682.000       +14092.000       +15213.000       +20227.000       +15203.000       +19005.000 ]      |            +46205.482
]

Hadamard比率9.2297e-30

delta=0.3 格基约化 =
[
        [          -1.000           +1.000           +0.000           +2.000           +2.000           +4.000           +4.000           +6.000 ]      |                +8.832
        [          +4.000           +4.000           +4.000           +3.000           +0.000           -1.000           -2.000           -3.000 ]      |                +8.426
        [          -1.000           -3.000           -3.000           +1.000           -1.000           -2.000           +1.000           +0.000 ]      |                +5.099
        [          -1.000           +1.000           +2.000           +1.000           +2.000           -2.000           +4.000           +1.000 ]      |                +5.657
        [          +4.000           +2.000           +0.000           -1.000           +2.000           -2.000           +0.000           -1.000 ]      |                +5.477
        [          +2.000           -2.000           +2.000           +1.000           -2.000           +2.000           +0.000           +2.000 ]      |                +5.000
        [          -2.000           +0.000           +0.000           +0.000           -3.000           -2.000           +1.000           +4.000 ]      |                +5.831
        [          +1.000           +0.000           +0.000           +1.000           -5.000           -1.000           -1.000           -3.000 ]      |                +6.164
]

Hadamard比率0.137705

正交性：
[ +1.00000  -0.32250  -0.13323  +0.44035  -0.28941  +0.31704  +0.31069  -0.64288 ]
[ -0.32250  +1.00000  -0.58187  +0.04196  +0.56336  +0.07121  -0.40706  +0.36579 ]
[ -0.13323  -0.58187  +1.00000  -0.03467  -0.32225  -0.11767  +0.33634  +0.19089 ]
[ +0.44035  +0.04196  -0.03467  +1.00000  +0.12910  -0.17678  +0.24254  -0.43015 ]
[ -0.28941  +0.56336  -0.32225  +0.12910  +1.00000  -0.25560  -0.43836  -0.05923 ]
[ +0.31704  +0.07121  -0.11767  -0.17678  -0.25560  +1.00000  +0.20580  +0.16222 ]
[ +0.31069  -0.40706  +0.33634  +0.24254  -0.43836  +0.20580  +1.00000  +0.05564 ]
[ -0.64288  +0.36579  +0.19089  -0.43015  -0.05923  +0.16222  +0.05564  +1.00000 ]

delta=0.75 格基约化 =
[
        [          +1.000           -1.000           -2.000           -1.000           -2.000           +2.000           -4.000           -1.000 ]      |                +5.657
        [          +1.000           +3.000           +3.000           -1.000           +1.000           +2.000           -1.000           +0.000 ]      |                +5.099
        [          +2.000           -2.000           +2.000           +1.000           -2.000           +2.000           +0.000           +2.000 ]      |                +5.000
        [          -4.000           -2.000           +0.000           +1.000           -2.000           +2.000           +0.000           +1.000 ]      |                +5.477
        [          -2.000           +0.000           +0.000           +0.000           -3.000           -2.000           +1.000           +4.000 ]      |                +5.831
        [          -2.000           -1.000           +1.000           +4.000           +2.000           +0.000           +0.000           +1.000 ]      |                +5.196
        [          +1.000           +0.000           +0.000           +1.000           -5.000           -1.000           -1.000           -3.000 ]      |                +6.164
        [          -2.000           +3.000           -2.000           +2.000           -1.000           +1.000           +3.000           +1.000 ]      |                +5.745
]

Hadamard比率0.343309

正交性：
[ +1.00000  -0.03467  +0.17678  +0.12910  -0.24254  -0.40825  +0.43015  -0.36927 ]
[ -0.03467  +1.00000  +0.11767  -0.32225  -0.33634  -0.15097  -0.19089  -0.10242 ]
[ +0.17678  +0.11767  +1.00000  +0.25560  +0.20580  +0.07698  +0.16222  -0.20889 ]
[ +0.12910  -0.32225  +0.25560  +1.00000  +0.43836  +0.38650  +0.05923  +0.28604 ]
[ -0.24254  -0.33634  +0.20580  +0.43836  +1.00000  +0.06601  +0.05564  +0.35825 ]
[ -0.40825  -0.15097  +0.07698  +0.38650  +0.06601  +1.00000  -0.34341  +0.20101 ]
[ +0.43015  -0.19089  +0.16222  +0.05923  +0.05564  -0.34341  +1.00000  -0.05648 ]
[ -0.36927  -0.10242  -0.20889  +0.28604  +0.35825  +0.20101  -0.05648  +1.00000 ]

delta=0.95 格基约化 =
[
        [          -1.000           -3.000           -3.000           +1.000           -1.000           -2.000           +1.000           +0.000 ]      |                +5.099
        [          +2.000           -2.000           +2.000           +1.000           -2.000           +2.000           +0.000           +2.000 ]      |                +5.000
        [          +4.000           +2.000           +0.000           -1.000           +2.000           -2.000           +0.000           -1.000 ]      |                +5.477
        [          -2.000           +0.000           +0.000           +0.000           -3.000           -2.000           +1.000           +4.000 ]      |                +5.831
        [          +1.000           +0.000           -1.000           +2.000           -1.000           -2.000           -3.000           +3.000 ]      |                +5.385
        [          +2.000           +1.000           -1.000           -4.000           -2.000           +0.000           +0.000           -1.000 ]      |                +5.196
        [          -2.000           +3.000           -2.000           +2.000           -1.000           +1.000           +3.000           +1.000 ]      |                +5.745
        [          -1.000           -1.000           +1.000           +5.000           -3.000           -1.000           -1.000           -2.000 ]      |                +6.557
]

Hadamard比率0.339015

正交性：
[ +1.00000  -0.11767  -0.32225  +0.33634  +0.21851  -0.15097  +0.10242  +0.29907 ]
[ -0.11767  +1.00000  -0.25560  +0.20580  +0.22283  -0.07698  -0.20889  +0.21350 ]
[ -0.32225  -0.25560  +1.00000  -0.43836  +0.03390  +0.38650  -0.28604  -0.36195 ]
[ +0.33634  +0.20580  -0.43836  +1.00000  +0.44585  -0.06601  +0.35825  +0.10461 ]
[ +0.21851  +0.22283  +0.03390  +0.44585  +1.00000  -0.21442  -0.09698  +0.28318 ]
[ -0.15097  -0.07698  +0.38650  -0.06601  -0.21442  +1.00000  -0.20101  -0.46957 ]
[ +0.10242  -0.20889  -0.28604  +0.35825  -0.09698  -0.20101  +1.00000  +0.10619 ]
[ +0.29907  +0.21350  -0.36195  +0.10461  +0.28318  -0.46957  +0.10619  +1.00000 ]