/*
  Author:  Pate Williams (c) 1997

  "Algorithm 2.6.3 (LLL Algorithm). Given
  a basis b[1], b[2],..., b[n] of a lattice
  (L, q), this algorithm transforms the vectors
  b[i] so that when the algorithm terminates the
  b[i] form an LLL-reduced basis." -Henri Cohen-
  See "A Course in Computational Algebraic
  Number Theory" by Henri Cohen pages 87-88.
  Solution of the subset problem using LLL
  reduction. See "Handbook of Applied
  Cryptography" by Alfred J. Menezes et al
  pages 120 - 121.
*/

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>

void system_error(char error_message[])
{
  fprintf(stderr, "%s",error_message);
  exit(1);
}

float **allocate_real_matrix(int lr, int ur, int lc, int uc)
{
  /* Allocates a real matrix of range [lr..ur][lc..uc]. */
  float **p = calloc((ur - lr + 1), sizeof(float *));
  int i;

  if (!p) system_error("Failure in allocate_real_matrix().");
  p -= lr;
  for (i = lr; i <= ur; i++) {
    p[i]= calloc((uc - lc + 1), sizeof(float));
	 if (!p[i]) system_error("Failure in allocate_real_matrix().");
    p[i] -= lc;
  }
  return p;
}

float *allocate_real_vector(int l, int u)
{
  /* Allocates a real vector of range [l..u]. */
  float *p;

  p = calloc((u - l + 1), sizeof(float));
  if (!p) system_error("Failure in allocate_real_vector().");
  return p - l;
}

void free_real_matrix(float **m, int lr, int ur, int lc)
{
  /* Frees a real matrix of range [lr..ur][lc..uc]. */
  int i;

  for (i = ur; i >= lr; i--) free((char *)(m[i] + lc));
  free((char *) (m + lr));
}

void free_real_vector(float *v, int l)
{
  /* Frees a real vector of range [l..u]. */

  free((char *)(v + l));
}

float Scalar(long n, float *u, float *v)
{
  /* Calculate the scalar product of two vectors
     [1..n] */
  double sum = 0.0;
  int i;

  for (i = 1; i <= n; i++) sum += u[i] * v[i];
  return sum;
}

void RED(int k, int l, int n,
         float **b, float **h,
         float **mu)
{
  int i, q = 0.5 + mu[k][l];

  if (fabs(mu[k][l]) > 0.5) {
    for (i = 1; i <= n; i++) {
      b[k][i] -= q * b[l][i];
      h[i][k] -= q * h[i][l];
    }
    mu[k][l] -= q;
    for (i = 1; i <= l - 1; i++) mu[k][i] -= q * mu[l][i];
  }
}

void SWAP(int k, int k1, int kmax, int n,
          float m, float *B, float **b,
          float **bs, float **h, float **mu)
{
  float C, t;
  float *c = allocate_real_vector(1, n);
  int i, j;

  for (j = 1; j <= n; j++) {
    t = b[k][j], b[k][j] = b[k1][j], b[k1][j] = t;
    t = h[j][k], h[j][k] = h[j][k1], h[j][k1] = t;
  }
  if (k > 2)
    for (j = 1; j <= k - 2; j++)
      t = mu[k][j], mu[k][j] = mu[k1][j], mu[k1][j] = t;
  C = B[k] + m * m * B[k1];
  mu[k][k1] = m * B[k1] / C;
  for (i = 1; i <= n; i++) c[i] = bs[k1][i];
  for (i = 1; i <= n; i++)
    bs[k1][i] = bs[k][i] + m * c[i];
  for (i = 1; i <= n; i++)
    bs[k][i] = - m * bs[k][i] + B[k] * c[i] / C;
  B[k] *= B[k1] / C;
  B[k1] = C;
  for (i = k + 1; i <= kmax; i++) {
    t = mu[i][k];
    mu[i][k] = mu[i][k1] - m * t;
    mu[i][k1] = t + mu[k][k1] * mu[i][k];
  }
  free_real_vector(c, 1);
}

void LLL(int n, float **b, float **h)
{
  /* Lattice reduction algorithm. */
  float m;
  float *B = allocate_real_vector(1, n);
  float **bs = allocate_real_matrix(1, n, 1, n);
  float **mu = allocate_real_matrix(1, n, 1, n);
  int i, j, k, k1, kmax = 1, l;

  for (i = 1; i <= n; i++) bs[1][i] = b[1][i];
  B[1] = Scalar(n, bs[1], bs[1]);
  for (i = 1; i <= n; i++) {
    for (j = 1; j <= n; j++)
      h[i][j] = 0.0;
    h[i][i] = 1.0;
  }
  for (k = 2; k <= n; k++) {
    if (k <= kmax) goto L3;
    kmax = k;
    for (i = 1; i <= n; i++) bs[k][i] = b[k][i];
    for (j = 1; j <= k - 1; j++) {
      mu[k][j] = Scalar(n, b[k], bs[j]) / B[j];
      for (i = 1; i <= n; i++)
        bs[k][i] -= mu[k][j] * bs[j][i];
    }
    B[k] = Scalar(n, bs[k], bs[k]);
    if (B[k] == 0.0)
      system_error("Failure in LLL.");
    L3:
    k1 = k - 1;
    m = mu[k][k1];
    RED(k, k1, n, b, h, mu);
    if (B[k] < (0.75 - m * m) * B[k1]) {
      SWAP(k, k1, kmax, n, m, B, b, bs, h, mu);
      k = max(2, k1);
      goto L3;
    }
    for (l = k - 2; l >= 1; l--)
      RED(k, l, n, b, h, mu);
  }
  free_real_matrix(bs, 1, n, 1);
  free_real_matrix(mu, 1, n, 1);
  free_real_vector(B, 1);
}

int SubsetSum(int n, float s, float *a, float *x)
{
  int n1 = n + 1;
  float **b = allocate_real_matrix(1, n1, 1, n1);
  float **h = allocate_real_matrix(1, n1, 1, n1);
  float sum;
  int i, j, m = ceil(sqrt(n) / 2.0);

  for (i = 1; i <= n1; i++) {
    if (i <= n) {
      b[i][i] = 1.0;
      b[i][n1] = m * a[i];
    }
    else {
      for (j = 1; j <= n; j++) b[i][j] = 0.5;
      b[i][n1] = m * s;
    }
  }
  printf("the matrix to be reduced is:\n\n");
  for (i = 1; i <= n1; i++) {
    for (j = 1; j <= n1; j++)
      printf("%6.2f ", b[i][j]);
    printf("\n");
  }
  printf("\n");
  LLL(n1, b, h);
  printf("the reduced matrix is:\n\n");
  for (i = 1; i <= n1; i++) {
    for (j = 1; j <= n1; j++)
      printf("%6.2f ", b[i][j]);
    printf("\n");
  }
  printf("\n");
  for (i = 1; i <= n1; i++) {
    for (j = 1; j <= n; j++) x[j] = b[i][j] + 0.5;
    sum = 0.0;
    for (j = 1; j <= n; j++) sum += a[j] * x[j];
    if (sum == s) {
      free_real_matrix(b, 1, n1, 1);
      free_real_matrix(h, 1, n1, 1);
      return 1;
    }
    for (j = 1; j <= n; j++) x[j] = - b[i][j] + 0.5;
    sum = 0.0;
    for (j = 1; j <= n; j++) sum += a[j] * x[j];
    if (sum == s) {
      free_real_matrix(b, 1, n1, 1);
      free_real_matrix(h, 1, n1, 1);
      return 1;
    }
  }
  free_real_matrix(b, 1, n1, 1);
  free_real_matrix(h, 1, n1, 1);
  return 0;
}

int main(void)
{
  int i, n = 8;
  float *a = allocate_real_vector(1, n);
  float *x = allocate_real_vector(1, n);
  float s;

  srand(time(NULL));
  printf("\n");
  for (i = 1; i <= n; i++) a[i] = pow(2, i);
  s = a[1 + rand() % n] + a[1 + rand() % n];
  if (SubsetSum(n, s, a, x)) {
    printf("sum: %f\n\n", s);
    printf("x[i]\t\ta[i]\n\n");
    for (i = 1; i <= n; i++)
      printf("%f\t%f\n", x[i], a[i]);
  }
  else printf("subset sum has no solution\n");
  free_real_vector(a, 1);
  free_real_vector(x, 1);
  return 0;
}