#if 0
Shamir Sharing
Warning!! We use the stock random number generator.
You must replace it if you really want to keep a secret!!

67 is prime and 67^4>2^24. We use the field of integers modulo 67.
We want to produce k<67 versions of a secret so that we may
reconstruct the message when q (0<q<k) versions of the secret are
available. ("q" for quorum)

We use polynomials of degree q-1 over the field.
If one knows values of p(x) for q distinct arguments
then the polynomial may be determined.

To code the secret we let each 6 bit byte of the message have its own
polynomial of degree q-1, but we speak here as if the secret were
just one character long. The secret character is s. We choose the
polynomial p such that p(0) = s. We choose the polynomial otherwise
to have random (mod 67) coefficients.
For 0<j<67, version j of the secret is p(j).

Suppose for some set S of q integers we know p(i) for each i in S.
To compute s = p(0) we apply Lagrange's interpolation formula:

p(x) = sum (i in S) (p(i) *
    (product (j in S but not i) (x - j))/
    (product (j in S but not i) (i - j)))

All arithmetic is modulo 67.
We use the special case where x=0.

The routines below deal in secrets < 67^4 which conveniently codes
three arbitrary bytes. Each resulting version is four bytes but
each of those bytes is < 67. Adding '0' to each version byte
results in portable ascii characters.

Program logic: "u25" is a type that gives 25 bit (or more)
unsigned integers.
Routine "void it(void)" initializes multiplication and division tables
modulo 67. mt[i][j] is i*j mod 67 if 0<=i<67 and 0<=j<67.
j*dt[j] = 1 mod 67 if  0<j<67.

Routine "void split(u25 sec, int q, int d, char w[d][4])"
splits the secret (sec) into d parts. A quorum of size q
is necessary to recover the secret. The parts are placed in
the array w. w[0] will hold part 1, and w[quor-1] will hold
part quor. As the parts are distributed to the custodians,
each part number must be included!

Routine "u25 join(int q, char w[q][4], int c[q])"
takes quor parts and recomputes the secret which it returns
as a function value. w[0] thru w[q-1] are the parts
and c[0] thru c[q-1] are the respective part numbers.

Routine main provides a fairly general test invocation of
this code.

The following stuff is missing from <stdlib.h> and libraries
in some systems claiming to be ANSI C!

static ldiv_t ldiv(long a, long b)
{ldiv_t A; A.quot = a/b; A.rem = a % b; return A;}
#endif

#define N 67
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
typedef unsigned long u25 /* 32 bits */;
typedef unsigned char uc;
static char mt[N][N], dt[N];
static void it(void){int i, j;
  for(i=0; i<N; ++i) for(j=0; j<N; ++j) {uc t = i*j % N; mt[i][j] = t;
      if(t==1) dt[j] = i;}}

static char M(char a){if(a<0) return a+N; return a;}

static u25 join(int quor, char w[quor][4], int c[quor])
{char sx[4]; int k; for(k=0; k<4; ++k) {char q = 0;
   {int n;   for(n=0; n<quor; ++n) {char p = w[n][k];
       {int m; for(m=0; m<quor; ++m)
         if(c[n]-c[m]) p = mt[p][mt[N-c[m]][dt[M(c[n]-c[m])]]];}
       q = M(q + p - N);}}
   sx[k]=q;}
return sx[3] + N*(sx[2] + N*(sx[1] + (u25)N*sx[0]));}

static void split(u25 sec, int quor, int dis, char w[dis][4])
{if(sec >= (u25)N*N*N*N) printf("Foul value\n");
 if(quor > dis || dis >= N)
   printf("Committee size must not exceed distribution and "
      "distribution must be less than N\n");
 {ldiv_t A = ldiv(sec, N), B = ldiv(A.quot, N), C = ldiv(B.quot, N);
  char a[4]; a[3] = A.rem; a[2] = B.rem; a[1] = C.rem; a[0] = C.quot;
  {int k; for(k = 0; k<4; ++k)
   {char coef[quor]; coef[0] = a[k]; 
     {int m; for(m = 1; m<quor; ++m) coef[m] = rand() % N;}
     {int n; for(n = 1; n<=dis; ++n) {char q = 0, m;
       for(m=quor-1; m>=0; --m) q = M(coef[m] + mt[q][n] - N);
       w[n-1][k] = q;}}}}}}

void sj(int d, int Q, int q[Q])
{char dx[d][4]; split((u25)12345678, Q, d, dx);
 printf("The %d shares; quorum = %d:\n", d, Q);
 {int i, j; for(i=0; i<d; ++i) {printf("%2d ", i+1);
   for (j=0; j<4; ++j) printf("%c", '0' + dx[i][j]); printf("\n");}}
 {char dz[Q][4]; {int i, j;
   for(i=0; i<Q; ++i) for(j=0; j<4; ++j) dz[i][j] = dx[q[i]-1][j];}
   printf("%ld\n", join(Q, dz, q));}}

int main(){it();
{int i, j; for(i=1; i<N; ++i) for(j=0; j<N; ++j)
   if(mt[i][mt[j][dt[i]]] != j) printf("Ouch, %d, %d\n", i, j);}
{int q[] = {1}; sj(2, 1, q);}
{int q[] = {2, 3}; sj(3, 2, q);}
{int q[] = {8, 3, 5, 1}; sj(8, 4, q);}
{int q[] = {8, 3, 4, 1}; sj(8, 4, q);}
{int q[] = {8, 3, 5, 3}; sj(8, 4, q);} // fails for duplicate portions.
{int q[] = {8, 3, 5, 1, 6}; sj(8, 5, q);}
{int q[] = {8, 4, 5, 1}; sj(9, 4, q);}
return 0;}