Blog literacki, portal erotyczny - seks i humor nie z tej ziemi
/*
Author: Pate Williams (c) 1997
Solution of the subset problem using LLL
reduction. See "Handbook of Applied
Cryptography" by Alfred J. Menezes et al
pages 120 - 121.
*/
#include
#include
#include
#include
void system_error(char error_message[])
{
printf("%s",error_message);
exit(1);
}
double **allocate_real_matrix(long m, long n)
{
long i;
double **p = calloc(m, sizeof(double *));
if (!p) system_error("Failure in allocate_real_matrix().");
for (i = 0; i < m; i++){
p[i]= calloc(n, sizeof(double));
if (!p[i]) system_error("Failure in allocate_real_matrix().");
}
return p;
}
double *allocate_real_vector(long n)
{
double *p = calloc(n, sizeof(double));
if (!p) system_error("Failure in allocate_real_vector().");
return p;
}
void free_real_matrix(double **m, long n)
{
long i;
for (i = 0; i < n; i++) free(m[i]);
free(m);
}
void free_real_vector(double *v)
{
free(v);
}
double Scalar(long n, double *u, double *v)
{
double sum = 0.0;
long i;
for (i = 0; i < n; i++) sum += u[i] * v[i];
return sum;
}
void Reduce(long k, long l, long n, double **b, double **mu)
{
long i, j, r = 0.5 + mu[k][l];
if (fabs(mu[k][l]) > 0.5) {
for (i = 0; i < n; i++) b[k][i] -= r * b[l][i];
for (j = 0; j < l; j++) mu[k][j] -= r * mu[l][j];
mu[k][l] -= r;
}
}
int LLL(long n, double **b)
{
/* Lattice reduction algorithm. */
double *B = allocate_real_vector(n);
double **bs = allocate_real_matrix(n, n);
double **mu = allocate_real_matrix(n, n);
double C, t, temp, x, y;
long i, j, k, l;
for (i = 0; i < n; i++) bs[0][i] = b[0][i];
B[0] = Scalar(n, bs[0], bs[0]);
for (i = 1; i < n; i++) {
for (j = 0; j < n; j++) bs[i][j] = b[i][j];
for (j = 0; j < i; j++) {
mu[i][j] = Scalar(n, b[i], bs[j]) / B[j];
for (k = 0; k < n; k++)
bs[i][k] -= mu[i][j] * bs[j][k];
}
B[i] = Scalar(n, bs[i], bs[i]);
}
L3:
k = 1;
L4:
l = k - 1;
Reduce(k, l, n, b, mu);
x = mu[k][l];
y = 0.75 - x * x;
if (B[k] < y * B[l]) {
C = B[k] + x * x * B[l];
mu[k][l] = x * B[l] / C;
B[k] *= B[l] / C;
B[l] = C;
for (i = 0; i < n; i++) {
temp = b[k][i];
b[k][i] = b[l][i];
b[l][i] = temp;
}
if (k > 1) {
for (j = 0; j < k - 1; j++) {
temp = mu[k][j];
mu[k][j] = mu[l][j];
mu[l][j] = temp;
}
}
for (i = k + 1; i < n; i++) {
t = mu[i][k];
mu[i][k] = mu[i][l] - x * t;
mu[i][l] = t + mu[k][l] * mu[i][k];
}
k = max(1, k - 1);
goto L4;
}
for (l = k - 2; l >= 0; l--) Reduce(k, l, n, b, mu);
k++;
if (k < n) goto L4;
free_real_matrix(bs, n);
free_real_matrix(mu, n);
free_real_vector(B);
return 1;
}
int SubsetSum(long n, double s, double *a, double *x)
{
long n1 = n + 1;
double **b = allocate_real_matrix(n1, n1);
double sum;
long i, j, m = ceil(sqrt(n) / 2.0);
for (i = 0; i < n1; i++) {
if (i < n) {
for (j = 0; j < n1; j++) b[i][j] = 0.0;
b[i][i] = 1.0;
b[i][n1 - 1] = m * a[i];
}
else {
for (j = 0; j < n; j++) b[i][j] = 0.5;
b[i][n1 - 1] = m * s;
}
}
printf("the matrix to be reduced is:\n\n");
for (i = 0; i < n1; i++) {
for (j = 0; j < n1; j++)
printf("%6.2f ", b[i][j]);
printf("\n");
}
printf("\n");
if (!LLL(n1, b)) {
free_real_matrix(b, n1);
return 0;
}
printf("the reduced matrix is:\n\n");
for (i = 0; i < n1; i++) {
for (j = 0; j < n1; j++)
printf("%6.2f ", b[i][j]);
printf("\n");
}
printf("\n");
for (i = 0; i < n1; i++) {
for (j = 0; j < n; j++) x[j] = b[i][j] + 0.5;
sum = 0.0;
for (j = 0; j < n; j++) sum += a[j] * x[j];
if (sum == s) {
free_real_matrix(b, n1);
return 1;
}
for (j = 0; j < n; j++) x[j] = - b[i][j] + 0.5;
sum = 0.0;
for (j = 0; j < n; j++) sum += a[j] * x[j];
if (sum == s) {
free_real_matrix(b, n1);
return 1;
}
}
free_real_matrix(b, n1);
return 0;
}
int main(void)
{
long i, n = 8;
double *a = allocate_real_vector(n);
double *x = allocate_real_vector(n);
double s;
srand(time(NULL));
printf("\n");
for (i = 0; i < n; i++) a[i] = pow(2, i);
s = a[rand() % n] + a[rand() % n];
if (SubsetSum(n, s, a, x)) {
printf("sum: %f\n\n", s);
printf("x[i]\t\ta[i]\n\n");
for (i = 0; i < n; i++)
printf("%f\t%f\n", x[i], a[i]);
}
else printf("subset sum has no solution\n");
free_real_vector(a);
free_real_vector(x);
return 0;
}