/* Functions to manipulate mod-2 polynomials */ #include #include #include /* 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: * 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; ndegree; 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; }