废话不多说，大名鼎鼎的Lenstra-Lenstra-Lovasz（LLL) 算法。实现参考论文：Factoring Polynomials with Rational Coefficients, 作者 A.K. Lenstra, H.W. Lenstra, Jr. and L. Lovasz.
/* File : LLL Lattice Reduction Matlab mex interface
* Description : do Lattice Reduction in the Real Field, Complex Lattice Redution is not currently supported
* reference paper 'Factoring Ploynominals with Rational Coefficients ' A.K Lenstra, H. W. Lenstra, and L.Lovasz
* Author : fangying
* Date : 2014-01-12
* Modified :
*/ #include "stdio.h"
#include "math.h"
#include "stdlib.h" /*
* Description : this function is used to give the round interger of the input double variable
* Input : a double variable
* Output : a round interger of the given variable
*/
int roundinteger(double var)
{
return (var > ) ? (int)(var + 0.5) : (int)(var - 0.5);
} /*
* Description : this function is used to perform inner product of two given vector or array
* Input : pointer of vector 'b1' and 'b2' with their length 'm'
* Output : pointer of the output result
*/ double innerproduct(double *b1, double *b2, int m)
{
double sum = ; while (m--)
{
sum += b1[m] * b2[m];
}
return sum;
} /*
* Description : this function performs Gram-Schmidt on input Matrix
* Input : pointer of the input Matrix(in array format) 'bin',pointer of the output Matrix 'bout',
* the projection between inner-vectors,and dimension 'N*M'('N' for cloumn,'M' for row)
* Output : the GSO(Gram-Schmidt Othogonal) Matrix pointer
*/ void GramSchmidt(double *bin, double *bout, double *u, double *B, int M, int N)
{
int i, j, k;
double uTemp, projection; /* copy bin to bout */
for (j = ; j < M*N; j++)
{
bout[j] = bin[j];
} /* Euclid Square of the first column */
B[] = innerproduct(bout, bout, M); for (i = ; i < N; i++)
{
for (k = ; k < i; k++)
{
projection = innerproduct(bout + k*M, bin + i*M, M);
uTemp = projection / B[k];
u[(i - )*(N - ) + k] = uTemp; for (j = ; j < M; j++)
{
bout[i*M + j] -= bout[k*M + j] * uTemp;
}
}
/* calculate the next B[i]*/
B[i] = innerproduct(bout + i*M, bout + i*M, M);
}
} /*
* Description : this function is used the do the reduction procedure
* Input : the column vector 'b',
* Output : update the ...
*/
void reduction(double *b, double *u, int k, int l, int M, int N)
{
int j, i, r;
double uTemp, uAbs; uTemp = u[(k - )*(N - ) + l];
uAbs = fabs(*u);
//uAbs = (uTemp >= 0) ? uTemp : (-uTemp);
//uAbs = (uTemp > 0) ? uTemp : (-uTemp); if (uAbs > 0.5)
{
r = roundinteger(uTemp); /* update u(k,k-1) <= u(k,k-1) - r */
u[(k - )*(N - ) + l] -= r; /* update b(k) <= b(k) -r*b(k-1) */
for (i = ; i < M; i++)
{
b[k*M + i] -= r*b[l*M + i];
} /* for u(k,j) with j<k-1, update u(k,j) <= u(k,j) - r*u(k-1,j) */
for (j = ; j <= l - ; j++)
{
u[(k - )*(N - ) + j] -= r*u[(l - )*(N - ) + j];
}
}
} /*
* Description : this function is used the do the check procedure
* Input : 'B' the input column Euclid Squre, 'b' the input Matrix element
'u' the projection , 'k' the current subscript and dimension 'M*N'
* Output : update the ...
*/
void check(double *B, double *b, double *u, int k, int l, int M, int N)
{
int i, j;
double uTemp, Btemp;
double tmp; uTemp = u[(k - )*(N - ) + k - ]; /* this is u(k,k-1) */ /* c(k-1) <= b(k)
* c(k) <= b(k-1)
* c(i) <= b(i) for i != k,k-1
*/
Btemp = B[k] + uTemp * uTemp * B[k - ]; // Btemp = c*(k-1) /* update u(k,k-1) <= u(k,k-1)B[k-1]/C[k-1] ， this is v(k,k-1) */
u[(k - )*(N - ) + k - ] = uTemp*B[k - ] / Btemp; /* update B[k] <= C[k] */
B[k] = B[k - ] * B[k] / Btemp;
/* update B[k-1] <= C[k-1] */
B[k - ] = Btemp; /* exchange b[k] <=> b[k] */
for (j = ; j < M; j++)
{
tmp = b[k*M + j];
b[k*M + j] = b[(k - )*M + j];
b[(k - )*M + j] = tmp;
} /* for j<k-1 ,i.e j = [1 to k-2]
* u(k-1,j) <= u(k,j)
* u(k,j) <= u(k-1,j)
*/
for (j = ; j < k - ; j++)
{
tmp = u[(k - )*(N - ) + j]; // u(k-1,j)
u[(k - )*(N - ) + j] = u[(k - )*(N - ) + j];
u[(k - )*(N - ) + j] = tmp; // u(k,j)
} /* for j>k
*
* u(i,k) <= u(i,k-1) - u(i,k)*u(k,k-1)
* u(i,k-1) <= u(k,k-1) -uTemp*u(i,k)
*/ for (i = k + ; i < N; i++)
{
tmp = u[(i - )*(N - ) + k]; // u(i,k)
/* v(i,k) <= u(i,k-1) - u(i,k)*u(k,k-1) */
u[(i - )*(N - ) + k] = u[(i - )*(N - ) + (k - )] - u[(i - )*(N - ) + k] * uTemp;
/* v(i,k-1) <= u(i,k-1)*v() */ //更新v(i,k-1)，看文献Page521页
u[(i - )*(N - ) + k - ] = u[(k - )*(N - ) + (k - )] * u[(i - )*(N - ) + (k - )] + tmp*( - uTemp*u[(k - )*(N - ) + (k - )]);
} } /*
* Description : LLL (Lattice Reduction Alogorithm
* Input : the input Matrix in array 'bin',the dimension of the Matrix 'N*M',the output 'HLR'
*/
void RLLL(double *bin, int M, int N, double *HLR)
{
int i, k, l; double *u;
double *B;
double *H; u = (double *)calloc(N*M, sizeof(double));
B = (double *)calloc(N, sizeof(double));
H = (double *)calloc(N*M, sizeof(double)); for (i = ; i < M*N; i++)
{
H[i] = bin[i];
} GramSchmidt(H, HLR, u, B, M, N); k = ; while ()
{
l = k - ;
reduction(H, u, k, l, M, N); /* iteration procedure */
if (B[k]<(0.75 - u[(k - )*(N - ) + k - ] * u[(k - )*(N - ) + k - ])*B[k - ])
{
check(B, H, u, k, l, M, N);
if (k>)
{
k--;
}
}
else
{
for (l = k - ; l >= ; l--)
{
reduction(H, u, k, l, M, N);
} if (k == N - ) break;
k++;
}
} for (i = ; i < M*N; i++)
{
HLR[i] = H[i];
} free(u);
free(H);
free(B); } int main()
{
int col = ;
int row = ;
int index = ; double magic[] = { , , , , , , , , };
double hlr[ * ] = { }; LLL(magic, col, row, hlr); return ;
}