Lattice Reduction（LLL）算法在C语言中实现

Lattice Reduction（LLL）算法在C语言中实现

采用了仿照手算的方法编写，使用了 “取出向量——>向量计算——>放回矩阵” 的模式进行计算，欢迎大家交流指正！

参考文献：Factoring Polynomials with Rational Coefficients, 作者 A.K. Lenstra, H.W. Lenstra, Jr. and L. Lovasz.

#include "stdafx.h"
#include <math.h>
#include "stdlib.h"
#include "stdio.h"
double delta = 0.75;//lovasz常数
double intermulti(double *A, double *B,int n)//正交计算
{
double sum = 0;
for (int i = 0; i < n; i++)
{
sum += A[i] * B[i];
}
return sum;
}
void add(double *IN1, double *IN2, double *OUT, int n)//向量加法
{
for (int i = 0; i < n; i++)
{
OUT[i] = IN1[i] + IN2[i];
}
}
void minus(double *IN1, double *IN2, double *OUT, int n)//向量减法
{
for (int i = 0; i < n; i++)
{
OUT[i] = IN1[i] - IN2[i];
}
}
void multiply(double IN1, double *IN2, double *OUT, int n)//向量点乘
{
for (int i = 0; i < n; i++)
{
OUT[i] = IN1*IN2[i];
}
}
double mui(double *A, double *B, int n)//求<a,b>/||b||^2过程
{
double a, b;
a = intermulti(A, B, n);
b = intermulti(B, B, n);
return a / b;
}
void vetor(int num, double *IN, double *OUT, int col,int row)//取出向量操作
{															//num为取出矩阵中向量序号
if (num > row) return;									
for (int i = 0; i <col; i++)
{
OUT[i] = IN[(num - 1)*col + i];
}
}
void Vetor_to_Array(double *Ve, double *Ar, int col, int row, int n)//将放回计算好的向量放回矩阵
{																	//n为将向量放在第n个位置
for (int i = 0; i < col; i++)									
{																					
Ar[(n - 1)*col + i] = Ve[i];
}
}
void Schmidt(double *IN, double *OUT, int col, int row)//施密特正交化
{
double k, *a, *b; int q;
a = (double *)calloc(col, sizeof(double));
b = (double *)calloc(col, sizeof(double));
for (int i = 0; i < col*row; i++)
{
if (i/col)
{
k = IN[i];
vetor(i / col + 1, IN, a, col, row);
for (q = 0; q < i/col; q++)
{
vetor(q + 1, OUT, b, col, row);
k -= mui(a,b,col)*OUT[q+(i%col)];
}
OUT[i] = k;
}
else
{
OUT[i] = IN[i];
}
}
free(a);
free(b);
}
void swap_vetor(double *IN, int col, int row, int A, int B)//矩阵中两向量交换位置
{															//A、B为两向量编号
double temp;
for (int i = 0; i < col; i++)
{
temp = IN[(A-1)*col + i];
IN[(A - 1)*col + i] = IN[(B - 1)*col + i];
IN[(B - 1)*col + i] = temp;
}
}
int Is_Fail_lovaszCondition(double *IN,double *IN_sch, int col, int row,int K)//判断是否符合lovasz条件，K无用可舍去
{												
double *a, *b, *c, *d; int count = 0;
a = (double *)calloc(col, sizeof(double));
b = (double *)calloc(col, sizeof(double));
c = (double *)calloc(col, sizeof(double));
d = (double *)calloc(col, sizeof(double));
for (int i = 1; i < row; i++)
{
vetor(i + 1, IN, a, col, row);
vetor(i, IN_sch, b, col, row);
vetor(i + 1, IN_sch, c, col, row);
multiply(mui(a, b, col), b, d, col);
add(d, c, d, col);
if (delta*intermulti(b, b, col)>intermulti(d, d, col))//lovasz条件判断
{
++count;
swap_vetor(IN, col, row, i, i + 1);
}
}
free(a);
free(b);
free(c);
free(d);
return count;
}
void Vetor_Change(double *Ar, double *Ar_sch, int I1, int J1, int col, int row)//Size-Reduction变换
{												     
double *a, *b,*b_sch, *d;
double c;
a = (double *)calloc(col, sizeof(double));
b = (double *)calloc(col, sizeof(double));
b_sch = (double *)calloc(col, sizeof(double));
d = (double *)calloc(col, sizeof(double));
vetor(I1, Ar, a, col, row);
vetor(J1, Ar, b, col, row);
vetor(J1, Ar_sch, b_sch, col, row);
c = round(mui(a,b_sch,col));
multiply(c, b, d, col);
minus(a, d, d, col);
Vetor_to_Array(d, Ar, col, row, I1);
free(a);
free(b);
free(b_sch);
free(d);
}
void LLL(double *IN, double *OUT, int col, int row)//LLL算法本体
{
double *a, *a_sch;
int k = 1, i, by = 0, qjr = 1;
a = (double *)calloc(col*row, sizeof(double));
a_sch = (double *)calloc(col*row, sizeof(double));
for ( i = 0; i <col*row ; i++)
{
a[i] = IN[i];
}
while (by+qjr)
{
Schmidt(a, a_sch, col, row);
for ( i = 1; i <= row; i++)
{
for (int j = i - 1; j >= 0; j--)
{
Vetor_Change(a, a_sch, i, j, col, row);
}
}
if (!Is_Fail_lovaszCondition(a, a_sch, col, row, k))
{
//k = k - 1 > 1 ? k - 1 : 1;
break;
}
}
for ( i = 0; i < col*row; i++)
{
OUT[i] = a[i];
}
free(a);
free(a_sch);
}
int _tmain(int argc, _TCHAR* argv[])//main函数
{
int col =  ,//列数
row =  ;//行数
double input[/*填入总元素数*/] = { /*填入你的矩阵，行向量形式*/ };
double output[/*填入总元素数*/] = { 0 };
LLL(input, output, col, row);
for (int i = 0; i < col*row; i++)
{
printf("%f ", output[i]);
if(i%col==(col-1)) printf("\n");
}
return 0;
}

本文地址：https://blog.csdn.net/qq_33127317/article/details/107217766