#include <gmp.h>
#include <md5.h>
#include <stdio.h>
#include <errno.h>
#include <stdarg.h>
#include <stdlib.h>
#include <unistd.h>
#include <strings.h>
#include <sys/wait.h>
#include <sys/time.h>
#include <sys/resource.h>

extern const char *__progname;

#define MAX_PARTS 4096

static int sflag = 0;
static int mflag = 0;
static int iflag = 0;
static int Iflag = 0;
static int obase = 10;
static int vflag = 0;
#define VLEVEL(n) (vflag >= (n))

static int argc;
static char **argv;

static int rbits;
static int rpool_n;
static unsigned char *rpool;
static int rpool_p;
#define RPOOL_UNREVEALED 256

static void printv(int, const char *, ...)
	__attribute__ ((__format__(__printf__,2,3)));
static void printv(int vlevel, const char *fmt, ...)
{
 va_list ap;

 if (! VLEVEL(vlevel)) return;
 va_start(ap,fmt);
 vfprintf(stderr,fmt,ap);
 va_end(ap);
}

static void handleargs(int ac, char **av)
{
 int errs;
 long int v;
 char *ep;

 errs = 0;
 argc = 0;
 argv = malloc((ac+1)*sizeof(char *));
 for (ac--,av++;ac;ac--,av++)
  { if (**av == '-')
     { for (++*av;**av;++*av)
	{ switch (**av)
	   { default:
		fprintf(stderr,"%s: bad flag -%c\n",__progname,**av);
		errs ++;
		break;
	     case 'm':
		mflag = 1;
		sflag = 0;
		break;
	     case 's':
		sflag = 1;
		mflag = 0;
		break;
	     case 'i':
		iflag = 1;
		break;
	     case 'I':
		Iflag = 1;
		break;
	     case 'v':
		vflag ++;
		break;
	     case 'x':
		obase = 16;
		break;
	     case 'b':
		++ *av;
		v = strtol(*av,&ep,0);
		if (*ep || (ep == *av))
		 { fprintf(stderr,"%s: invalid number %s\n",__progname,*av);
		   for (ep=*av;*ep;ep++) ;
		   errs ++;
		 }
		else
		 { obase = v;
		 }
		*av = ep - 1;
		break;
	   }
	}
     }
    else
     { argv[argc++] = *av;
     }
  }
 if (errs) exit(1);
}

static void checkargs(void)
{
 int errs;

 errs = 0;
 if (mflag)
  { if ((obase < 2) || (obase > 36))
     { fprintf(stderr,"%s: output base must be [2..36]\n",__progname);
       errs ++;
     }
    if (argc < 2)
     { fprintf(stderr,"%s: need at least two parts with -m\n",__progname);
       errs ++;
     }
  }
 else if (sflag)
  { if (argc != 4)
     { fprintf(stderr,"%s: need exactly four args with -s\n",__progname);
       errs ++;
     }
  }
 else
  { fprintf(stderr,"%s: need either -s or -m\n",__progname);
    errs ++;
  }
 if (errs) exit(1);
}

static void random_stir(void)
{
 void *mv;
 int o;
 unsigned char r[16];
 int i;

 for (o=(rpool_n-1U)&~15U;o>=0;o-=16)
  { mv = md5_init();
    md5_process_bytes(mv,rpool,rpool_n);
    md5_result(mv,&r[0]);
    for (i=0;i<16;i++) rpool[o+i] ^= r[i];
  }
}

static void randomness_from_data(const void *buf, int nb)
{
 const unsigned char *bp;

 bp = buf;
 for (;nb>0;nb--)
  { rpool[rpool_p++] ^= *bp++;
    if (rpool_p >= rpool_n)
     { random_stir();
       rpool_p = RPOOL_UNREVEALED / 8;
     }
  }
}

static void randomness_from_fd(int fd)
{
 char buf[8192];
 int nb;

 while (1)
  { nb = read(fd,&buf[0],sizeof(buf));
    randomness_from_data(&nb,sizeof(nb));
    if (nb <= 0) break;
    randomness_from_data(&buf[0],nb);
  }
}

static void randomness_from_command(const char *execfile, ...)
{
 int nargs;
 const char **args;
 const char *arg;
 va_list ap;
 int p[2];
 int pid;
 const char *env[1];

 va_start(ap,execfile);
 for (nargs=0;va_arg(ap,const char *);nargs++) ;
 va_end(ap);
 args = malloc((nargs+1)*sizeof(const char *));
 va_start(ap,execfile);
 for (nargs=0;(arg=va_arg(ap,const char *));nargs++) args[nargs] = arg;
 va_end(ap);
 args[nargs] = 0;
 env[0] = 0;
 if (pipe(p) < 0)
  { fprintf(stderr,"%s: pipe: %s\n",__progname,strerror(errno));
    exit(1);
  }
 pid = fork();
 if (pid < 0)
  { fprintf(stderr,"%s: fork: %s\n",__progname,strerror(errno));
    exit(1);
  }
 if (pid == 0)
  { if (p[1] != 1)
     { dup2(p[1],1);
       close(p[1]);
     }
    close(p[0]);
    execve(execfile,(const void *)args,(const void *)env);
    _exit(1);
  }
 close(p[1]);
 randomness_from_fd(p[0]);
 close(p[0]);
 while ((pid = wait3(0,WNOHANG,0)) > 0) ;
}

static void quick_randomness(void)
{
 { struct timeval tv;
    gettimeofday(&tv,0);
    randomness_from_data(&tv,sizeof(tv));
  }
  { struct rusage ru;
    getrusage(RUSAGE_SELF,&ru);
    randomness_from_data(&ru,sizeof(ru));
    getrusage(RUSAGE_CHILDREN,&ru);
    randomness_from_data(&ru,sizeof(ru));
  }
  { unsigned long int i;
    i = getpid();
    randomness_from_data(&i,sizeof(i));
    i = getppid();
    randomness_from_data(&i,sizeof(i));
    i = getuid();
    randomness_from_data(&i,sizeof(i));
  }
}

static void init_random(void)
{
 printv(1,"Collecting randomness");
 rpool_n = ((rbits + 127U + RPOOL_UNREVEALED) & ~127U) / 8U;
 rpool = malloc(rpool_n);
 rpool_p = 0;
 /* deliberately don't init it */
 if (! Iflag)
  { printv(2,"\n (ps awwlxe)");
    randomness_from_command("/bin/ps","ps","awwlxe",(char *)0);
    printv(2,"\n (w)");
    randomness_from_command("/usr/bin/w","w",(char *)0);
    printv(2,"\n (ls -ali /tmp/.)");
    randomness_from_command("/bin/ls","ls","-ali","/tmp/.",(char *)0);
    printv(2,"\n (netstat -s)");
    randomness_from_command("/usr/bin/netstat","netstat","-s",(char *)0);
    printv(2,"\n (netstat -na)");
    randomness_from_command("/usr/bin/netstat","netstat","-na",(char *)0);
    printv(2,"\n (netstat -ni)");
    randomness_from_command("/usr/bin/netstat","netstat","-ni",(char *)0);
  }
 if (iflag || Iflag)
  { printv(2,"\n (stdin)");
    randomness_from_fd(0);
  }
 printv(2,"\n (quick sources)");
 quick_randomness();
 printv(2,"\n (stirring)");
 random_stir();
 printv(1,"\n");
}

static unsigned char random_byte(void)
{
 unsigned char rv;

 rv = rpool[rpool_p++];
 if (rpool_p >= rpool_n)
  { random_stir();
    rpool_p = RPOOL_UNREVEALED / 8;
  }
 return(rv);
}

static void random_integer(MP_INT *num, unsigned int nbits)
{
 unsigned int bytes;
 char *strform;
 int i;

 bytes = (nbits + 7) / 8;
 strform = malloc((bytes*2)+1);
 for (i=0;i<bytes;i++) sprintf(strform+i+i,"%02x",random_byte());
 mpz_set_str(num,strform,16);
 bzero(strform,bytes*2);
 mpz_mod_2exp(num,num,nbits);
}

static void get_arglist_number(MP_INT *val, const char *s0)
{
 int base;
 char *ep;
 const char *s;

 base = 0;
 if (s0[0] == ':')
  { base = strtol(s0,&ep,10);
    if ((ep == s) || (base < 2) || (base > 36))
     { fprintf(stderr,"%s: invalid base in number %s\n",__progname,s);
       exit(1);
     }
    s = (*ep == ':') ? ep+1 : s0;
  }
 else
  { s = s0;
  }
 if (mpz_set_str(val,s,base) < 0)
  { fprintf(stderr,"%s: invalid number %s\n",__progname,s0);
    exit(1);
  }
}

static void compose_ints(MP_INT *rv, int n, ...)
{
 int b;
 int bb;
 int bit;
 int j;
 int v;
 int any;
 va_list ap;

 mpz_set_ui(rv,0);
 bb = 0;
 b = 0;
 do
  { bit = 1 << b;
    any = 0;
    va_start(ap,n);
    for (j=0;j<n;j++)
     { v = va_arg(ap,int);
       if (v >= bit) any = 1;
       if (v & bit) mpz_setbit(rv,bb);
       bb ++;
     }
    va_end(ap);
    b ++;
  } while (any);
}

static void decompose_ints(MP_INT *iv, int n, ...)
{
 int bit;
 int i;
 int vn;
 int *vv;
 va_list ap;
 MP_INT t;

 vv = malloc(n*sizeof(int));
 bzero(vv,n*sizeof(int));
 mpz_init(&t);
 mpz_set(&t,iv);
 bit = 1;
 vn = 0;
 while (mpz_sgn(&t) > 0)
  { if (mpz_get_ui(&t) & 1) vv[vn] |= bit;
    mpz_tdiv_q_2exp(&t,&t,1);
    vn ++;
    if (vn >= n)
     { vn = 0;
       bit <<= 1;
     }
  }
 va_start(ap,n);
 for (i=0;i<n;i++) *va_arg(ap,int *) = vv[i];
 va_end(ap);
}

static void split_it(void)
{
 MP_INT secret;
 MP_INT modulus;
 int totalparts;
 int minparts;
 int mbits;
 MP_INT aux;
 MP_INT *coeffs;
 int i;
 int p;

 mpz_init(&aux);
 get_arglist_number(&aux,argv[2]);
 if ((mpz_cmp_ui(&aux,2) <= 0) || (mpz_cmp_ui(&aux,MAX_PARTS) > 0))
  { fprintf(stderr,"%s: totalparts out of range [2..%d]\n",__progname,MAX_PARTS);
    exit(1);
  }
 totalparts = mpz_get_ui(&aux);
 get_arglist_number(&aux,argv[3]);
 if ((mpz_cmp_ui(&aux,2) <= 0) || (mpz_cmp_ui(&aux,totalparts) > 0))
  { fprintf(stderr,"%s: minparts out of range [2..%d]\n",__progname,totalparts);
    exit(1);
  }
 minparts = mpz_get_ui(&aux);
 mpz_init(&secret);
 mpz_init(&modulus);
 get_arglist_number(&secret,argv[0]);
 get_arglist_number(&modulus,argv[1]);
 if (mpz_cmp(&secret,&modulus) >= 0)
  { fprintf(stderr,"%s: secret must be less than modulus\n",__progname);
    exit(1);
  }
 mbits = mpz_sizeinbase(&modulus,2);
 rbits = (minparts+1) * mbits;
 init_random();
 printv(1,"Choosing polynomial coefficients\n");
 coeffs = malloc(minparts*sizeof(MP_INT));
 for (i=0;i<minparts-1;i++)
  { mpz_init(coeffs+i);
    random_integer(&aux,mbits+64);
    mpz_mod(coeffs+i,&aux,&modulus);
    if (VLEVEL(2))
     { fprintf(stderr,"  coeffs[%d] = 0x",i);
       mpz_out_str(stderr,16,coeffs+i);
       fprintf(stderr,"\n");
     }
  }
 mpz_init(coeffs+i);
 mpz_set(coeffs+i,&secret);
 if (VLEVEL(2))
  { fprintf(stderr,"  coeffs[%d] = 0x",i);
    mpz_out_str(stderr,16,coeffs+i);
    fprintf(stderr,"\n");
  }
 for (p=1;p<=totalparts;p++)
  { printf("1z");
    mpz_out_str(stdout,35,&modulus);
    printf("z");
    compose_ints(&aux,3,p,minparts,totalparts);
    mpz_out_str(stdout,35,&aux);
    printf("z");
    mpz_set_ui(&aux,0);
    for (i=0;i<minparts;i++)
     { mpz_mul_ui(&aux,&aux,p);
       mpz_add(&aux,&aux,coeffs+i);
       mpz_mod(&aux,&aux,&modulus);
       if (VLEVEL(3))
	{ fprintf(stderr,"p=%d i=%d value = 0x",p,i);
	  mpz_out_str(stderr,16,&aux);
	  fprintf(stderr,"\n");
	}
     }
    if (VLEVEL(2) && !VLEVEL(3))
     { fprintf(stderr,"poly(%d) = ",p);
       mpz_out_str(stderr,16,&aux);
       fprintf(stderr,"\n");
     }
    mpz_out_str(stdout,35,&aux);
    printf("\n");
  }
}

static void invalid_part(int, const char *, ...)
	__attribute__ ((__format__(__printf__,2,3)));
static void invalid_part(int pno, const char *whyfmt, ...)
{
 va_list ap;

 fprintf(stderr,"%s: invalid part %s: ",__progname,argv[pno]);
 va_start(ap,whyfmt);
 vfprintf(stderr,whyfmt,ap);
 va_end(ap);
 fprintf(stderr,"\n");
 exit(1);
}

static void parts_mismatch(const char *what)
{
 fprintf(stderr,"%s: parts don't agree with one another: %s mismatch\n",__progname,what);
 exit(1);
}

static void panic(const char *, ...)
	__attribute__ ((__format__(__printf__,1,2)));
static void panic(const char *fmt, ...)
{
 va_list ap;

 fprintf(stderr,"%s: THIS CAN'T HAPPEN: ",__progname);
 va_start(ap,fmt);
 vfprintf(stderr,fmt,ap);
 va_end(ap);
 fprintf(stderr,"\n");
 abort();
}

static void dumpmatrix(MP_INT **m, MP_INT *b, int nr, int nc)
{
 int r;
 int c;
 int w;
 int *cw;

 cw = malloc((nc+1)*sizeof(int));
 for (c=0;c<=nc;c++) cw[c] = 0;
 for (r=0;r<nr;r++)
  { for (c=0;c<nc;c++)
     { w = mpz_sizeinbase(&m[r][c],16);
       if (w > cw[c]) cw[c] = w;
     }
    w = mpz_sizeinbase(&b[r],16);
    if (w > cw[nc]) cw[nc] = w;
  }
 for (r=0;r<nr;r++)
  { for (c=0;c<nc;c++)
     { w = mpz_sizeinbase(&m[r][c],16);
       fprintf(stderr,"%*s",cw[c]+1-w,"");
       mpz_out_str(stderr,16,&m[r][c]);
     }
    w = mpz_sizeinbase(&b[r],16);
    fprintf(stderr," |%*s",cw[nc]+1-w,"");
    mpz_out_str(stderr,16,&b[r]);
    fprintf(stderr,"\n");
  }
 free(cw);
}

static void merge_it(void)
{
 MP_INT modulus;
 int totalparts;
 int minparts;
 int mbits;
 int tp;
 int mp;
 int pp;
 MP_INT aux1;
 MP_INT aux2;
 int p;
 char *sp;
 char *cp;
 int *pnos;
 MP_INT *b;
 MP_INT **m;
 MP_INT *m0;
 int i;
 int c;
 int r;

 mpz_init(&modulus);
 mpz_init(&aux1);
 mpz_init(&aux2);
 for (p=0;p<argc;p++)
  { if ((argv[p][0] != '1') || (argv[p][1] != 'z')) invalid_part(p,"unrecognized format value");
    sp = argv[p] + 2;
    cp = index(sp,'z');
    if (! cp) invalid_part(p,"missing modulus terminator");
    *cp = '\0';
    if ((mpz_set_str(&aux2,sp,35) < 0) || (mpz_sgn(&aux2) <= 0))
     { *cp = 'z';
       invalid_part(p,"invalid modulus");
     }
    *cp = 'z';
    sp = cp + 1;
    cp = index(sp,'z');
    if (! cp) invalid_part(p,"missing small values terminator");
    *cp = '\0';
    if ((mpz_set_str(&aux1,sp,35) < 0) || (mpz_sgn(&aux1) <= 0))
     { *cp = 'z';
       invalid_part(p,"invalid small values");
     }
    *cp = 'z';
    sp = cp + 1;
    decompose_ints(&aux1,3,&pp,&mp,&tp);
    if (VLEVEL(3))
     { fprintf(stderr,"small values = ");
       mpz_out_str(stderr,16,&aux1);
       fprintf(stderr," = pp=%d mp=%d tp=%d\n",pp,mp,tp);
     }
    if ((mpz_set_str(&aux1,sp,35) < 0) || (mpz_sgn(&aux1) < 0))
     { invalid_part(p,"invalid part value");
     }
    if (mpz_cmp(&aux1,&aux2) >= 0) invalid_part(p,"part value out of range");
    if (p == 0)
     { if ((tp < 2) || (tp > MAX_PARTS)) invalid_part(p,"total part count out of range [2..%d]",MAX_PARTS);
       totalparts = tp;
       if ((mp < 2) || (mp > tp)) invalid_part(p,"minimum part count out of range [2..%d]",totalparts);
       minparts = mp;
       mpz_set(&modulus,&aux2);
       mbits = mpz_sizeinbase(&modulus,2);
       pnos = malloc(argc*sizeof(int));
       b = malloc(argc*sizeof(MP_INT));
       for (i=0;i<argc;i++) mpz_init(b+i);
     }
    else
     { if (tp != totalparts) parts_mismatch("total part count");
       if (mp != minparts) parts_mismatch("minimum part count");
       if (mpz_cmp(&aux2,&modulus) != 0) parts_mismatch("modulus");
     }
    if ((pp < 1) || (pp > totalparts)) invalid_part(p,"impossible part number");
    pnos[p] = pp;
    mpz_set(b+p,&aux1);
  }
 if (VLEVEL(1))
  { fprintf(stderr,"modulus = ");
    mpz_out_str(stderr,16,&modulus);
    fprintf(stderr,"\n");
  }
 for (p=1;p<argc;p++)
  { for (i=p-1;i>=0;i--)
     { if (pnos[p] == pnos[i])
	{ if (mpz_cmp(b+p,b+i) == 0)
	   { fprintf(stderr,"%s: part %s repeated\n",__progname,argv[i]);
	   }
	  else
	   { fprintf(stderr,"%s: parts %s and %s are inconsistent\n",__progname,argv[i],argv[p]);
	   }
	  exit(1);
	}
     }
  }
 if (argc < minparts)
  { fprintf(stderr,"%s: too few parts supplied, can't solve\n",__progname);
    exit(1);
  }
 m = malloc(argc*sizeof(MP_INT *));
 m0 = malloc(argc*minparts*sizeof(MP_INT));
 m[0] = m0;
 for (i=(argc*minparts)-1;i>=0;i--) mpz_init(&m[0][i]);
 for (i=1;i<argc;i++) m[i] = m[i-1] + minparts;
 for (p=0;p<argc;p++)
  { mpz_set_ui(&aux1,1);
    for (i=minparts-1;i>0;i--)
     { mpz_set(&m[p][i],&aux1);
       mpz_mul_ui(&aux1,&aux1,pnos[p]);
       mpz_mod(&aux1,&aux1,&modulus);
     }
    mpz_set(&m[p][0],&aux1);
  }
 printv(1,"Parts are superficially consistent, matrix set up.\n");
 if (VLEVEL(2)) dumpmatrix(m,b,argc,minparts);
 /* At this point we need to solve M X = B (M is a matrix, X a column vector
    of the unknowns, and B a column vector of the knowns).  Since we're in a
    finite field instead of the reals, any nonzero value will do equally well
    as a pivot.  If we can't find _any_ pivot, something is wrong - this
    "can't happen", since we checked that the pnos[] values are all different
    and thus all the equations are known to be linearly inependent.  We
    checked that argc >= minparts, and thus M has at least as many rows as
    columns.  If more, we check that the extras are redundant... */
 for (c=0;c<minparts;c++)
  { if (VLEVEL(3))
     { fprintf(stderr,"In loop, matrix now\n");
       dumpmatrix(m,b,argc,minparts);
     }
    for (i=c;i<argc;i++) if (mpz_sgn(&m[i][c]) != 0) break;
    if (i >= argc) panic("matrix is singular");
    if (i != c)
     { MP_INT *tp;
       MP_INT tmi;
       printv(3,"Swapping rows %d and %d to obtain nonzero pivot\n",i,c);
       tp = m[i];
       m[i] = m[c];
       m[c] = tp;
       tmi = b[i];
       b[i] = b[c];
       b[c] = tmi;
     }
    mpz_gcdext(&aux2,&aux1,0,&m[c][c],&modulus);
    if (mpz_cmp_ui(&aux2,1) != 0)
     { fprintf(stderr,"%s: invalid split - modulus isn't prime: gcd(",__progname);
       mpz_out_str(stderr,10,&modulus);
       fprintf(stderr,",");
       mpz_out_str(stderr,10,&m[c][c]);
       fprintf(stderr,") = ");
       mpz_out_str(stderr,10,&aux2);
       fprintf(stderr,"\n");
       exit(1);
     }
    mpz_mod(&aux1,&aux1,&modulus);
    if (VLEVEL(3))
     { fprintf(stderr,"Pivot value is 0x");
       mpz_out_str(stderr,16,&m[c][c]);
       fprintf(stderr,"\nReciprocal is 0x");
       mpz_out_str(stderr,16,&aux1);
       fprintf(stderr,"\n");
     }
    for (i=0;i<c;i++) if (mpz_sgn(&m[c][i]) != 0) panic("nonzero value in row");
    for (i=c;i<minparts;i++)
     { mpz_mul(&m[c][i],&m[c][i],&aux1);
       mpz_mod(&m[c][i],&m[c][i],&modulus);
     }
    if (mpz_cmp_ui(&m[c][c],1) != 0) panic("reciprocal wasn't");
    mpz_mul(&b[c],&b[c],&aux1);
    mpz_mod(&b[c],&b[c],&modulus);
    for (r=0;r<argc;r++)
     { if (r == c) continue;
       for (i=c+1;i<minparts;i++)
	{ mpz_mul(&aux1,&m[c][i],&m[r][c]);
	  mpz_sub(&aux1,&m[r][i],&aux1);
	  mpz_mod(&m[r][i],&aux1,&modulus);
	}
       mpz_mul(&aux1,&b[c],&m[r][c]);
       mpz_sub(&aux1,&b[r],&aux1);
       mpz_mod(&b[r],&aux1,&modulus);
       mpz_set_ui(&m[r][c],0);
     }
  }
 if (VLEVEL(2))
  { fprintf(stderr,"Elimination done, matrix now\n");
    dumpmatrix(m,b,argc,minparts);
  }
 for (r=minparts;r<argc;r++)
  { if (mpz_sgn(&b[r]) != 0)
     { fprintf(stderr,"%s: invalid split - redundant parts disagree\n",__progname);
       exit(1);
     }
  }
 mpz_out_str(stdout,obase,&b[minparts-1]);
 printf("\n");
}

int main(int, char **);
int main(int ac, char **av)
{
 handleargs(ac,av);
 checkargs();
 if (sflag) split_it();
 if (mflag) merge_it();
 exit(0);
}