/*
  Author:  Pate Williams (c) 1997

  Exercise
  "5.10 Suppose the Merkle-Hellman Cryptosystem has
  as its public list of sizes the vector
  t = (1394, 1256, 1508, 1987, 439, 650, 339,
       2303, 810).
  Suppose Oscar discovers that p = 2503.
  (a) By trial and error, determine the value a
      such that a ^ -1 t mod p is a permutation
      of a superincreasing list.
  (b) Show how the ciphertext 5746 would be
      decrypted." -Douglas R. Stinson-
  See "Cryptography: Theory and Practice" by
  Douglas R. Stinson.
*/

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

long Extended_Euclidean(long b, long n)
{
  long b0 = b, n0 = n, t = 1, t0 = 0, temp, q, r;

  q = n0 / b0;
  r = n0 - q * b0;
  while (r > 0) {
    temp = t0 - q * t;
    if (temp >= 0) temp = temp % n;
    else temp = n - (- temp % n);
    t0 = t;
    t = temp;
    n0 = b0;
    b0 = r;
    q = n0 / b0;
    r = n0 - q * b0;
  }
  if (b0 != 1) return 0;
  else return t % n;
}

int superincreasing(long n, long *s)
{
  long i, j, sum;

  for (j = 0; j < n; j++) {
    sum = 0;
    for (i = 0; i < j; i++) sum += s[i];
    if (s[j] <= sum) return 0;
  }
  return 1;
}

int subset_sum(long T, long n, long *s, long *x)
{
  long S = T, i, sum = 0;

  for (i = n - 1; i >= 0; i--) {
    if (S >= s[i]) {
      S -= s[i];
      x[i] = 1;
    }
    else x[i] = 0;
  }
  for (i = 0; i < n; i++)
    sum += x[i] * s[i];
  return sum == T;
}

int main(void)
{
  int found = 0;
  long a, af, ai, i, j, n = 10, p = 2503, u, s[10];
  long x[10], y = 5746, z;
  long t[10] = {1394, 1256, 1508, 1987, 439,
                650, 724, 339, 2303, 810};

  for (a = 1; !found && a < p; a++) {
    ai = Extended_Euclidean(a, p);
    for (i = 0; i < n; i++)
      s[i] = (ai * t[i]) % p;
    /* sort the sizes into ascending order */
    for (i = 0; i < n - 1; i++)
      for (j = i + 1; j < n; j++)
         if (s[i] > s[j])
           u = s[i], s[i] = s[j], s[j] = u;
    found = superincreasing(n, s);
    af = a;
  }
  printf("the public list of sizes is:\n");
  for (i = 0; i < n; i++)
    printf("%4ld ", t[i]);
  printf("\n");
  printf("a = %ld a ^ - 1 = %ld\n", af, ai);
  printf("the superincreasing size vector is:\n");
  for (i = 0; i < n; i++)
    printf("%4ld ", s[i]);
  printf("\n");
  z = (ai * y) % p;
  printf("z = (%ld * %ld) %% %ld = %ld\n",
         ai, y, p, z);
  if (subset_sum(z, n, s, x)) {
    printf("the solution vector is:\n");
    for (i = 0; i < n; i++)
      printf("%ld ", x[i]);
    printf("\n");
  }
  else printf("subset sum has no solution\n");
  return 0;
}