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

Consider the finite field GF(256).
We want to produce k(<256) 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.
The values of p(x) for q distinct x's determine the polynomial.

To code the secret we let each 8 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 256) coefficients.
For 0<j<256, 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 in GF(256).
We use the special case where x=0.

Program logic: "ui" is a type that gives 32 bits.
Routine "void it(void)" initializes multiplication and division tables
for GF(256). mt[i][j] (or Mt(i, j)) is i*j.
j*dt[j] = 1 if  0<j<256.

Routine "void split(ui sec, int q, int d, uc 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 "ui join(int q, uc 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!

#endif

#define md 0
// md is a performance option. clang and x86 prefer 0 today for speed.
// md=1 saves 64Kb of space.

#define N 256
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
typedef unsigned int ui /* 32 bits */;
typedef unsigned char uc;
static ui const irp = 27;
#if md
typedef struct{ui q; ui r;} PR;
ui pmul(ui a, ui b){ui sm=0;
  while(b) {if(b&1) sm ^= a; a = a<<1; b>>=1;} return sm;}
static int High(ui x){int y=0; // position of most significant bit.
{int k=4; while(k--) {int p = 1<<k; if(x&-1L<<p) {x >>= p; y += p;}}}
return y;}
static PR pqr(ui n, ui d){if(!d) exit(printf("foss\n")); // synthetic division
  int hn = High(n), hd = High(d); ui a=0;
  int k = hn-hd; while(k>=0) {if(n&(1L<<(hn--))) {n ^= d<<k; a |= 1L<<k;} --k;}
  return (PR){a, n};}
static PR ge(ui a, ui b) // Extended Euclidean Algorithm for GF(2^q)
// After http://cap-lore.com/code/Scheme/repository/egcd
{if (!b) return (PR){1, 0};
 PR qr = pqr(a, b);
 if(!qr.r) return (PR){0, 1};
 PR c = ge(b, qr.r);
 return (PR){c.r, c.q ^ pmul(c.r, qr.q)};}
static ui finv(ui x){return ge(256+irp, x).r;}
static uc mt[3*16*16], dt[N];
#else
static uc mt[N][N], dt[N];
#endif
// fmul is the field multiply.
static ui fmul(char a, ui b){ui ts=N/2, sm=0;
  while(b) {if(b&1) sm ^= a;
    a = (a<<1)^(1?((a>>31)&irp):(a&ts?irp:0)); b>>=1;}
  return sm&(N-1);}

#if md
static uc Mt(uc j, uc k){return 0? // 1 is performance option. 1 wins today.
   mt[(j&(256-16))|(k>>4)]
 ^ mt[256|(j&(256-16))|(k&15)]
 ^ mt[256|((j&15)<<4)|(k>>4)]
 ^ mt[512|((j&15)<<4)|(k&15)]
 : fmul(j, k);}
void it(void){for(int j=0; j<16; ++j) for(int k=0; k<16; ++k) {
  mt[(j<<4)|k] = fmul(j<<4, k<<4);
  mt[256+((j<<4)|k)] = fmul(j<<4, k);
  mt[512+((j<<4)|k)] = fmul(j, k);}
if(1) for(int j=0; j<16; ++j) for(int k=0; k<16; ++k) 
  if (fmul(j, k<<4) ^ mt[256+((j<<4)|k)]) exit(printf("failA %d %d\n", j, k));
for(int j=0; j<256; ++j) dt[j] = finv(j);
if(1) for(int j=1; j<256; ++j)
  if(Mt(j, dt[j]) ^ 1) exit(printf("failB %x\n", j));
if(1) for(int j=0; j<256; ++j) for(int k=0; k<256; ++k)
  if(fmul(j, k) ^ Mt(j, k)){
              exit(printf("FailC %d %d\n", j, k));}}
#else
static uc Mt(uc j, uc k){return mt[j][k];}
void it(void){for(int i=0; i<N; ++i) for(int j=0; j<N; ++j)
  {uc t = fmul(i, j); mt[i][j] = t; if(t==1) dt[j] = i;}}
#endif

ui join(int quor, uc w[quor][4], int c[quor])
{uc sx[4]; for(int k=0; k<4; ++k) {uc q = 0;
   for(int n=0; n<quor; ++n) {uc p = w[n][k];
     for(int m=0; m<quor; ++m)
       if(c[n]^c[m]) p = Mt(p, Mt(c[m], dt[c[n]^c[m]]));
     q = q ^ p;}
   sx[k]=q;}
if(0) {for(int k=0; k<4; ++k) printf(":%02x", sx[k]); printf("\n");}
return sx[3] + N*(sx[2] + N*(sx[1] + (ui)N*sx[0]));}

void split(ui sec, int quor, int dis, uc w[dis][4])
{if(quor > dis || dis >= N)
   exit(printf("Quorum size (%d) must not exceed distribution (%d) and "
      "distribution must be less than %d\n", quor, dis, N));
 uc a[4]; {int j=4; while(j--) {a[j] = sec&255; sec >>= 8;}}
 for(int k = 0; k<4; ++k) {uc coef[quor]; coef[0] = a[k]; 
   for(int m = 1; m<quor; ++m) coef[m] = random() & 255;
   for(int n = 1; n<=dis; ++n) {uc q = 0;
     for(int m=quor-1; m>=0; --m) q = coef[m] ^ Mt(q, n);
     w[n-1][k] = q;}}}

static void sj(int d, int Q, int q[Q])
{ui sec = 123456789;
 uc dx[d][4]; split(sec, Q, d, dx);
 printf("%d shares; quorum = {", d);
   for(int j=0; j<Q; ++j) printf("%d ", q[j]); printf("}\n");
 for(int i=0; i<d; ++i) {printf("%2d ", i+1);
   for (int j=0; j<4; ++j) printf("%02x", dx[i][j]); printf("\n");}
 {uc dz[Q][4];
   for(int i=0; i<Q; ++i) for(int j=0; j<4; ++j) dz[i][j] = dx[q[i]-1][j];
   ui A = join(Q, dz, q); printf("%d %s\n\n", A, A == sec?"good":"bad");}}

#include <assert.h>
int main(){it();
  assert((char)-1 < 0); // lest fmul fail on machines with unsigned char's.
if(0){uc z = 190; // peform
  for(int j=0; j<8; ++j) z = Mt(z, z);
  printf("z = %x\n", z);
} else { // function
{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[] = {7, 3, 4, 1}; sj(8, 4, q);}
{int q[] = {8, 3, 5, 3}; sj(8, 4, q);} // fails for duplicate shares.
{int q[] = {8, 3, 5, 1, 6}; sj(8, 5, q);}
{int q[] = {8, 3, 5, 1}; sj(9, 4, q);}}
return 0;}