/* Functions to manipulate mod-2 polynomials  */

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


/* Data structures to store mod-2 polynomials
 *
 *  degree:	If the polynomial is nonzero, this is its nonnegative degree.
 *		If degree<0, the polynomial is the zero polynomial
 *
 *  coeff:	For the zero polynomial, this is NULL.  Otherwise, coeff[n]
 *		is the coefficient of t^n.  Note: coeff[n] is in {0,1}.
 *
 */
struct Mod2Poly {
  int degree;
  int *coeff;
};
typedef struct Mod2Poly *Mod2Poly;


/* **********************************************\
 * CREATE, ASSIGN, AND DESTROY MOD2POLY STRUCTS *
\************************************************/

/* ************************************************************
 * Allocate a Mod2Poly and its coeff[] member.  If coeffs==NULL
 * or degree<0, leave  new->coeff[]==NULL so the result is treated
 * as a zero polynomial. 
 */
Mod2Poly
newmod2poly( int *coeffs, int deg )
{
  Mod2Poly new;
  new = (Mod2Poly)malloc(sizeof(struct Mod2Poly)); assert(new);
  new->degree = deg;
  if( deg<0 ) {			/* ==> the zero polynomial */
    new->coeff = 0;
  } else {
    int n;
    new->coeff = (int *)calloc(1+deg, sizeof(int)); assert(new->coeff);
    if( coeffs )
      for(n=0; n<=deg; n++) new->coeff[n] = coeffs[n];
    /* ...else `new->coeff[]' will remain all zeroes. */
  }
  return new;
}

/* ************************************************************
 * Deallocate a Mod2Poly `p' and any non-NULL coeff[] member.
 * Check and do nothing if `p' is NULL.
 */
Mod2Poly
freemod2poly( Mod2Poly p )
{
  if(p) {
    if(p->coeff) free(p->coeff);
    free(p);
  }
  return 0;
}

/* ************************************************************
 * Assign the contents of one Mod2Poly into another.
 * If the target is NULL, allocate it.
 * If the target is too short, extend its coeff[] array.
 */
Mod2Poly
assignmod2poly( Mod2Poly target, Mod2Poly source )
{
  if(!target) {			/* ...then supply a new one */
    target = newmod2poly(source->coeff, source->degree);
  } else {			/* ...use the one provided */
    if ( target->degree != source->degree ) {
      if(source->degree<0) {	/* ==> `source' is the zero polynomial */
	if(target->coeff) free(target->coeff);
	target->coeff = 0;
	target->degree = -1;
      } else {
	int n;
	target->coeff = (int *)realloc(target->coeff, 
				       (1+source->degree)*sizeof(int));
	assert(target->coeff);
	for(n=0; n<=source->degree; n++) target->coeff[n] = source->coeff[n];
	target->degree = source->degree;
      }
    } else {			/* target->degree == source->degree */
      int n;
      for(n=0; n<=source->degree; n++) target->coeff[n] = source->coeff[n];
    }
  }
  return target;
}

/* ************************************************************
 * Restrict to the lowest-order terms of a Mod2Poly.
 * Compute the correct degree of the restriction.
 * Reallocate its `coeff[]' member if necessary.
 * Reassign its `degree' member, if necessary.
 * Update the Mod2Poly member values by reference.
 * Argument d0 controls the result as follows:
 *   d0 < 0                trim all of p to its actual degree
 *   d0 >= p->degree       trim all of p to its actual degree
 *   -1< d0 < p->degree    trim p->coeff[0],...,[d0] to its actual degree
 */
void
trimmod2poly( Mod2Poly p, int d0 )
{
  int d;
  if(d0<0 || d0>p->degree) d0=p->degree;
  if(p->degree>=0) {		/* If `p' is possibly nonzero... */
    assert(p->coeff);		/*   There must be a coefficients array */
    for(d = d0; d>=0 && p->coeff[d]==0; d--)
      ;				/*   Calculate the true degree. */
    p->degree = d;		/*   Set the true degree in the struct. */
    if(d>=0) {			/*   If true degree is nonnegative... */
      p->coeff = (int *)realloc(p->coeff, (1+d)*sizeof(int)); 
      assert(p->coeff);		/*     ...trim `coeff[]' in all cases. */
    } else {			/*   Else true degree is negative... */
      free(p->coeff);		/*     ...so eliminate `coeff[]'...  */
      p->coeff = 0;		/*     ...and leave a NULL member.  */
    }
  } else {			/* Else `p->degree' is negative... */
    if (p->coeff) {		/*   ...so make a zero polynomial... */
      free(p->coeff);		/*   ...with no coefficient array... */
      p->coeff = 0;		/*   ...and a null `coeff' member.  */
    }
  }
  return;
}

/* ************************************************************
 * Return the integer value of a Mod2Poly.
 */
unsigned int
mod2polytoint( Mod2Poly p )
{
  unsigned int val=0;
  if(p)
    if(p->degree>=0) {
      int n;
      for ( n=p->degree; n>=0; n-- )
	val = 2*val + p->coeff[n];
    }
  return val;
}

/* *****************************************************
 * Convert an integer to a Mod2Poly
 */
Mod2Poly
mod2polyfromint( unsigned int k )
{
  Mod2Poly out;
  int deg= -1, d;
  unsigned int k0;

  for(k0=k; k0>0; k0>>=1) ++deg; /* First find the degree `deg' */
  out = newmod2poly(0,deg);	/* Zero coefficients, to start */
  for( d=0; d<=deg; d++) {	/* Fill the coefficient array */
    out->coeff[d] = k&1;
    k>>=1;
  }
  return out;
}

/* *****************************************************
 * Print the contents of a Mod2Poly
 */
void
printmod2poly( FILE *outfile, Mod2Poly p, char *polyname )
{
  fprintf(outfile, "Mod-2 polynomial %s", polyname);
  if(!p) {
    fprintf(outfile, " is null.\n");
  } else {
    if(p->degree<0) {
      fprintf(outfile, " is the zero polynomial (degree p < 0).\n");
    } else {
      int n;
      fprintf(outfile, " has degree %d and coefficients:", p->degree);
      for (n=0; n<=p->degree; n++) {
	fprintf(outfile, " %d", p->coeff[n]);
      }
      fputc('\n', outfile);
    }
  }
  fflush(outfile);
}



/* ********************************************\
 * FUNCTIONS FROM MATHEMATICS FOR MULTIMEDIA. *
\* ********************************************/

/* ************************************************************
 * Mod-2 Polynomial Sums
 * 
 * mod2polysum( p[], dp, q[], dq ):
 * [0] Let ds = max(dp,dq) and allocate sum[0],...,sum[ds]
 * [1] For d=0 to ds, let sum[d] = 0
 * [2] For d=0 to dp, replace sum[d] ^= p[d]
 * [3] For d=0 to dq, replace sum[d] ^= q[d]
 * [4] Return sum[] and ds
 * 
 */
Mod2Poly
mod2polysum(Mod2Poly p, Mod2Poly q)
{
  Mod2Poly sum;
  int n;
  if(p->degree > q->degree) {
    sum = newmod2poly(p->coeff, p->degree);
    if(sum->coeff && q->coeff)
      for(n=0; n<=q->degree; n++) sum->coeff[n] ^= q->coeff[n];
  }
  else {
    sum = newmod2poly(q->coeff, q->degree);
    if(sum->coeff && p->coeff)
      for(n=0; n<=p->degree; n++) sum->coeff[n] ^= p->coeff[n];
  }
  trimmod2poly(sum, -1);
  return sum;
} 

/* ************************************************************
 * Mod-2 Polynomial Products
 * 
 * mod2polyproduct( p[], dp, q[], dq ):
 * [0] Let ds = dp+dq and allocate prod[0],...,prod[ds]
 * [1] For d=0 to ds, let prod[d] = 0
 * [2] For i=0 to dp, do [3]
 * [3]    For j=0 to dq, replace prod[i+j] ^= p[i]&q[j]
 * [4] Return prod[] and ds
 * 
 */
Mod2Poly
mod2polyproduct(Mod2Poly p, Mod2Poly q)
{
  Mod2Poly prod;
  int deg, i, j;

  if(p->degree<0 || q->degree<0) {
    prod = newmod2poly(0,-1);	/* Zero factor ==> zero product */
  }
  else {
    deg = p->degree + q->degree;
    prod = newmod2poly(0, deg);	/* All zero `coeff[]' to start */
    for(i=0; i<=p->degree; i++)
      for(j=0; j<=q->degree; j++)
	prod->coeff[i+j] ^= p->coeff[i] & q->coeff[j];
  }
  trimmod2poly(prod, -1);
  return prod;
} 

/* ************************************************************
 * Mod-2 Polynomial Degee
 * 
 * mod2polydegree( p[], dp ):
 * [0] For d=dp down to 0, do [1]
 * [1]   If p[d]==1, then return d
 * [2] Return -1
 * 
 * In this implementation, use d0 for the initial degree guess dp.
 *   d0 < 0                return actual degree of p
 *   d0 >= p->degree       return actual degree of p
 *   -1< d0 < p->degree    return degree of p->coeff[0],...,[d0]
 */
int
mod2polydegree(Mod2Poly p, int d0 )
{
  int d;
  if(d0<0 || d0>p->degree) d0=p->degree;
  for(d=d0; d>=0; d-- )
    if( p->coeff[d] ) return d;
  return -1;
} 

/* ************************************************************
 * Mod-2 Polynomial Quotient and Remainder
 * 
 * mod2polydivision( p1[], d1, p2[], d2 ):
 * [0] If d1<d2, then let dq = -1 and let dr = d1
 * [1] Else do [2] to [5]
 * [2]    Let dq = d1-d2 and let dr = d2-1
 * [3]    For d=d1 down to d2, do [4] to [5]
 * [4]       If p1[d]==1, then do [5]
 * [5]          For n=0 to d2-1, replace p1[n+d-d2] ^= p2[n]
 * [6] Let dr = mod2polydegree(p1[], dr)
 * [7] Return p1[], dq, dr
 * 
 * In this implementation, the function's return value is degree
 * of the remainder r, while  the coefficients of r and q are
 * returned by reference as parts of p1:
 *
 * If d1<d2:
 *   dr = d1
 *   dq = -1 and q=0
 *   p1[0],...,p1[dr] are the coefficients of r=p1
 *
 * If d1>=d2:
 *   dq = d1-d2
 *   If dr>=0
 *     r!=0
 *     p1[0],...,p1[dr] are the coefficients of r
 *     p1[d2],...,p1[d2+dq] are the coefficients of q
 *   Else dr<0
 *     r==0
 *     p1[0]=...=p1[d2-1]=0 are the coefficients of r=0
 *     p1[d2],...,p1[d2+dq] are the coefficients of q
 *
 */
int
mod2polydivision( Mod2Poly p1, Mod2Poly p2 )
{
  int dr, d, n;
  if(mod2polydegree(p2,-1)<0) 	/* exit; don't divide by 0 */
    return p1->degree;
  if(mod2polydegree(p2,-1)==0) 	/* exit; don't divide by 1 */
    return -1;
  if(p1->degree < p2->degree) {
    dr = p1->degree;
  } else {
    dr = p2->degree - 1;
    for(d=p1->degree; d>=p2->degree; d-- )
      if(p1->coeff[d])
	for(n=0; n<p2->degree; n++)
	  p1->coeff[n+d-p2->degree] ^= p2->coeff[n];
  }
  return mod2polydegree(p1, dr);
} 

/* ************************************************************
 * Euclid's Algorithm for Mod-2 Polynomials
 * 
 * mod2polygcd( x[], dx, y[], dy ):
 * [0-] Let dz = dx and allocate z[0],...,z[dz]
 * [0]  For d=0 to dx, copy z[d] = x[d]
 * [1]  Compute dr,dq,y[] with mod2polydivision(y[],dy,x[],dx)
 * [1+] Let dx = dr, and for d=0 to dr, copy x[d] = y[d]
 * [2]  Let dy = dz, and for d=0 to dz, copy y[d] = z[d]
 * [3]  If dr>=0, then go to [0]
 * [4]  Return y[], dy
 * 
 * Special cases like x==1 are handled separately.
 * Note that the quotient mod-2 polynomial is never used.
 */
Mod2Poly
mod2polygcd( Mod2Poly x, Mod2Poly y)
{
  Mod2Poly z=0;
  int dr=0;

  while( mod2polydegree(x,-1)>=0 ) {
    z = assignmod2poly( z, x );
    dr = mod2polydivision(y,x);
    if(dr<0) {			/* ...then x <- zero mod-2 polynomial */
      if(x->coeff) {
	free(x->coeff);
	x->coeff=0;
      }
      x->degree= -1;
    } else {			/* ...else x <- remainder part of y */
      trimmod2poly(y, dr);
      assignmod2poly( x, y );
    }
    assignmod2poly( y, z );	/* ...and y <- archived value of x */
  }
  freemod2poly(z);
  return y; 
}