/*****************************************************************************/
/*                                                                           */
/* Author: Michael Monagan, June 2024                                        */
/*                                                                           */
/* Code to factor a multilinear polynomial of F2                             */
/*                                                                           */
/* f = (abc+a+c+ab) (x+y+1)                                                  */
/*                                                                           */
/*****************************************************************************/
/*                                                                           */
/* MULMOD7.c   Added enlargepoly, mulpoly2, equalpoly, facpoly2              */
/* MULMOD8.c   Do two hash tests ==> Prob(error) <= n^2 deg f^2 / (2^63)^2   */
/* MULMOD10.c  Add setbit, resizepoly, and printing from 1, mergesort        */
/*                                                                           */
/*****************************************************************************/


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

#define LONG unsigned long int
#define ULONG unsigned long int

#define DEBUGfacpoly 0

ULONG * array64u( int n ) {
// allocate an array of n ULONG's not n bytes
    unsigned N;
    ULONG *A;
    if( n==0 ) return NULL;
    N = sizeof(ULONG) * n;
    A = (ULONG *) malloc(N);
    if( A==0 ) { printf("out of memory\n"); exit(1); }
    return A;
}

typedef struct poly {
    int n;  // number of variables
    int t;  // number of terms
    ULONG *terms; // bit vector
} poly;

/************************************************************************/
/*                                                                      */
/*  Encode a boolean polynomial in F2[x0,x1,x2] in an array of bits     */
/*  If n = 80 use two 64 bit words for each monomial thus 2 t words     */
/*  f = x0x1 + x0x2 + x1 + x2 is encoded as {3,4,[11,101,10,100]}       */
/*  g = x0x1x3 + x0x2 + x3 is encoded as {4,3,[1011,101,1000]}          */
/*                                                                      */
/************************************************************************/

int termsize(int n) { return (n+63)/64; }

void printpoly( poly *f, char x ) {
    int n,t,q,r,m,i,j,k,first; ULONG mask, *A, z;
    n = f->n;
    t = f->t;
    if( t==0 ) { printf("0;\n"); return; }
    A = f->terms;
    q = n/64;
    r = n%64;
    m = termsize(n);
    for( i=0; i<t; i++, A+=m ) {
        if( i>0 ) printf("+");
        first = 1;
        for( n=1,j=0; j<q; j++ ) { // printing from 1 not 0
            z = A[j];
            for( mask=1, k=0; k<64; k++, n++, mask<<=1 )
                if( z&mask ) { 
                    if(!first) printf("*");
                    printf("%c%d",x,n);
                    first = 0;
                }
        }
        if( r ) {
            z = A[j];
            for( mask=1, k=0; k<r; k++, n++, mask<<=1 )
                if( z&mask ) { 
                    if(!first) printf("*");
                    printf("%c%d",x,n);
                    first = 0;
                }
        }
        if( first ) printf("1");
    }
    printf(";\n");
    return;
}

void printpoly2( poly *f, char x ) {
// Allow polynomials to have degree < 4
    int n,t,q,r,m,i,j,k,first,d; ULONG mask, *A,z;
    n = f->n;
    t = f->t;
    if( t==0 ) { printf("0;\n"); return; }
    A = f->terms;
    q = n/64;
    r = n%64;
    m = termsize(n);
    for( i=0; i<t; i++,A+=m ) {
        if( i>0 ) printf("+");
        first = 1;
        for( n=0, j=0; j<q; j++ ) {
            z = A[j];
            for( mask=3, k=0; k<32; k++, n++, mask<<=2 )
                if( z&mask ) { 
                    if(!first) printf("*");
                    d = (z&mask) >> 2*k;
                    if( d==1 ) printf("%c%d",x,n); else printf("%c%d^%d",x,n,d);
                    first = 0;
                }
        }
        if( r ) {
            z = A[j];
            for( mask=3, k=0; k<r/2; k++, n++, mask<<=2 )
                if( z&mask ) { 
                    if(!first) printf("*");
                    d = (z&mask) >> 2*k;
                    if( d==1 ) printf("%c%d",x,n); else printf("%c%d^%d",x,n,d);
                    first = 0;
                }
        }
        if( first ) printf("1");
    }
    printf(";\n");
    return;
}


// constructor for poly
poly * newpoly(int n, int t) {
    int i,m;
    ULONG *A;
    poly *f;
    f = (poly *) malloc( sizeof(poly) );
    f->n = n;
    f->t = t;
    m = termsize(n);
    A = array64u(m*t);
    for( i=0; i<m*t; i++ ) A[i] = 0;
    f->terms = A;
    return f;
}


// destructor for poly
void freepoly(poly *f) {
    ULONG *A;
    A = f->terms;
    free(A);
    free(f);
}

// copy constructor
poly * copypoly( poly *f ) {
    int n,t,m,i,j; ULONG *A,*B;
    poly* g;
    n = f->n;
    t = f->t;
    A = f->terms;
    m = termsize(n);
    g = newpoly(n,t);
    B = g->terms;
    for( i=0; i<t; i++,A+=m,B+=m )
        for( j=0; j<m; j++ ) B[j] = A[j];
    return g;
}

poly * resizepoly(poly *f, int s) {
    int n,t,m,i; ULONG *A,*B; poly *g;
    n = f->n;
    t = f->t;
    A = f->terms;
    m = termsize(n);
    g = newpoly(n,s);
    B = g->terms;
    if( t>s ) t = s;
    for( i=0; i<m*t; i++ ) B[i] = A[i];
    freepoly(f);
    return g;
}


/*************************   listpoly   **********************************/

typedef struct listpoly {
    poly * f;
    struct listpoly *next;
} listpoly;


listpoly * conspoly(poly *f, listpoly * p) {
    listpoly * q;
    q = (listpoly *) malloc( sizeof(listpoly) );
    q->f = f;
    q->next = p;
    return q;
}

listpoly *joinpoly(listpoly *p, listpoly *q) {
    listpoly *r;
    if( p==NULL ) return q;
    if( q==NULL ) return p;
    r = p; while( r->next != NULL ) r = r->next;
    r->next = q;
    return p;
}


ULONG *indetspoly(poly *f) {
    int n,t,m,i,j;
    ULONG *A,*X;
    n = f->n;
    t = f->t;
    A = f->terms;
    m = termsize(n);
    X = array64u(m);
    for( j=0; j<m; j++ ) X[j] = 0;
    for( i=0; i<t; i++ ) {
        for( j=0; j<m; j++ ) X[j] = X[j] | A[j]; 
        A += m;
    }
    return X;
}

ULONG seed;
ULONG mult;

ULONG rand64u() {
// return 64 random bits
    ULONG x,y;
    extern ULONG seed, mult;
    seed = mult*seed;
    x = seed >> 32;
    seed = mult*seed;
    y = seed >> 32;
    x = (x<<32) | y;
    return(x);
}

void sortpoly(poly *f);
int compare(int m, ULONG *X, ULONG *Y);

ULONG * randbits( int m, int n ) {
// return X = array64u(n) where X[i] is random on [1,2^m) for 0<m<=64
    int i; ULONG mask,x,*X;
    if( m<1 || m>64 ) { printf("randbits: m must be in [1,64]\n"); exit(1); }
    if( n<1 ) { printf("randbits: n must be > 0\n"); exit(1); }
    X = array64u(n);
    mask = 0xFFFFFFFFFFFFFFFFL;
    mask = mask >> (64-m);
    for( i=0; i<n; i++ ) { 
        // do { x = rand64u() & mask; } while( x==0 );
        x = rand64u() & mask; X[i] = x;
    }
    return X;
}

poly* randpoly( int n, int t ) {
    int q,r,m,i,j,s;
    ULONG mask, *A;
    poly *f;
    f = newpoly(n,t);
    A = f->terms;
    q = n/64;
    r = n%64;
    if( r==0 ) m = q; else m = q+1;
    mask = 0xFFFFFFFFFFFFFFFFL;
    if( r ) mask = mask >> (64-r);
    for( i=0; i<t; i++ ) {
        for( j=0; j<q; j++ ) A[j] = rand64u();
        if( r ) A[j] = rand64u() & mask;
        A += m;
    }
    sortpoly(f);
    A = f->terms;
    for( i=0,s=0; i<t; i=j,s++ ) {
        if( s<i ) for( j=0; j<m; j++ ) A[s*m+j] = A[i*m+j];
        for( j=i+1; j<t && compare(m,A+m*i,A+m*j)==0; j++ ); 
    }
    f->t = s;
    return f;
}

poly* randpoly1( int n, int k, int t ) {
// set the first k bits of each monomial to random, zero out the remaining n-k bits
    int q,r,m,i,j,s;
    ULONG mask, *A;
    poly *f;
    m = termsize(n);
    f = newpoly(n,t);
    A = f->terms;
    for( i=0; i<t; i++,A+=m )
        for( j=0; j<m; j++ ) A[j] = 0;
    q = k/64;
    r = k%64;
    mask = 0xFFFFFFFFFFFFFFFFL;
    if( r ) mask = mask >> (64-r);
    A = f->terms;
    for( i=0; i<t; i++ ) {
        for( j=0; j<q; j++ ) A[j] = rand64u();
        if( r ) A[j] = rand64u() & mask;
        A += m;
    }
    sortpoly(f);
    A = f->terms;
    for( i=0,s=0; i<t; i=j,s++ ) {
        if( s<i ) for( j=0; j<m; j++ ) A[s*m+j] = A[i*m+j];
        for( j=i+1; j<t && compare(m,A+m*i,A+m*j)==0; j++ ); 
    }
    f->t = s;
    return f;
}

void shiftbits( ULONG *X, int m, int s ) {
    //  shift the bits in X right (up) by s>=0 bits
    int i,q; ULONG x,y,mask;
    if( s<0 ) { printf("shiftbits s must be nonnegative\n"); exit(1); }
    if( s>63 ) { 
        q = s/64;
        for( i=m-1; i>=q; i-- ) X[i] = X[i-q];
        while( i>=0 ) { X[i] = 0; i--; }
        s = s%64;
    }
    if( s==0 ) return;
    mask = 0xFFFFFFFFFFFFFFFFL;
    mask = mask << (64-s);
    for( i=m-1; i>=0; i-- ) {
        x = X[i] << s;
        if( i>0 ) { y = X[i-1] & mask; y = y >> (64-s); x = x | y; }
        X[i] = x;
    }
    return;
}
    
void shiftpoly( poly *f, int s ) {
    int n,t,m,i; ULONG *A;
    n = f->n;
    t = f->t;
    A = f->terms;
    m = termsize(n);
    for( i=0; i<t; i++,A+=m ) shiftbits(A,m,s);
    return;
}

/*
poly* randpoly2( int n, int k, int t ) {
// set the first k bits of each monomial to 0
    poly *f;
    f = randpoly(n,t);
    shiftpoly(f,k);
    return f;
}
*/

poly* randpoly2( int n, int k, int t ) {
// set the first k bits of each monomial to 0
    int q,r,m,i,j,s;
    ULONG mask, *A;
    poly *f;
    f = randpoly(n,t);
    m = termsize(n);
    q = k/64;
    r = k%64;
    mask = 0xFFFFFFFFFFFFFFFFL;
    if( r ) mask = ( mask >> r) << r;
    //printf("randpoly2: mask=%lu\n",mask);
    A = f->terms;
    for( i=0; i<t; i++,A+=m ) {
        for( j=0; j<q; j++ ) A[j] = 0;
        if( r ) A[j] = A[j] & mask;
    }
    sortpoly(f);
    A = f->terms;
    for( i=0,s=0; i<t; i=j,s++ ) {
        if( s<i ) for( j=0; j<m; j++ ) A[s*m+j] = A[i*m+j];
        for( j=i+1; j<t && compare(m,A+m*i,A+m*j)==0; j++ ); 
    }
    f->t = s;
    return f;
}

/************   HASHING HERE *************/
ULONG degpoly(ULONG m);
ULONG counter;

ULONG MULMOD( ULONG a, ULONG b, ULONG m, ULONG B, ULONG M ) {
// Given a,b,m in F_2[z] with deg(a)<d, deg(b)<d where d=deg(m)>0
// Compute a x b mod m in F_2[z] -- Rem(a*b,m,z) mod 2; in Maple
// B=2^deg(b) and M=2^d are bit masks
   ULONG c;
   counter++;
   if( a==0 || b==0 ) return 0;
   c = 0;
   while(B) {
      c = c << 1; // c = 2c
      if( c&M ) c = c ^ m; // c = c xor m
      if( b&B ) c = c ^ a; // c = c xor a
      B = B >> 1; // B = B/2
   }
   return c;
}

ULONG mulmodGF264( ULONG a, ULONG b, ULONG m, ULONG bbit ) {
// Let M = z^64 + a63 z^63 + ... + a1 z + a0 be in GF(2)[z]
// Given A,B in GF(2)[z] with deg(A)<64, deg(B)<64 compute AB mod M
// The A(z), B(z) and M(z) are input as a=A(2), B=B(2), m=M(2)-2^64
   ULONG c,topbit;
   counter++;
   if( a==0 || b==0 ) return 0;
   topbit = 1LU<<63;
   c = 0;
   while(bbit) {
      if( c&topbit ) c = (2*c)^m; else c = 2*c; // c = zc mod m
      if( b&bbit ) c = c^a; // c = c+a
      bbit = bbit >> 1;
   }
   return c;
}

ULONG mulmod( ULONG a, ULONG b, ULONG m) {
// Compute a x b mod m in Z2[z]  (Rem(a*b,m,z) mod 2)
    ULONG c,B,M;
    if( a==0 || b==0 ) return 0;
    if( a>=m ) { printf("a must be in [0,m)\n"); exit(1); }
    if( b>=m ) { printf("b must be in [0,m)\n"); exit(1); }
    B = 1L << degpoly(b);
    M = 1L << degpoly(m);
    c = MULMOD(a,b,m,B,M);
    return c;
}

ULONG hashprod(int n, ULONG *X, ULONG *H, ULONG minpoly, ULONG *B, ULONG M) {
// X is a bit vector of length n representing a monomial in n varialbes
// H is n (random) elements of F_2[z]/minpoly(z) encoded in an integer
// Output product( H[i]^X[i], i=0..n-1 ) mod minpoly(z)
    int q,r,i,j;
    ULONG h,mask;
    q = n/64;
    r = n%64;
    h = 1;
    for( i=0; i<q; i++ )
        for( mask=1, j=0; j<64; j++, mask<<=1 )
            if( X[i]&mask ) h = MULMOD( h, H[64*i+j], minpoly, B[64*i+j], M );
        for( mask=1, j=0; j<r; j++, mask<<=1 )
            if( X[i]&mask ) h = MULMOD( h, H[64*i+j], minpoly, B[64*i+j], M );
    return h;
}


ULONG degpoly(ULONG m) {
// if m = 101 corresponding to z^2+1 then return 2
    int d;
    for( d=0; m>1; d++ ) m = m>>1;
    return d;
}


int haspoly(poly *f,int x) {
// return true iff has(f,x) == degree(f,x)>0
    int n,t,m,q,r,i,j;
    ULONG *A,mask;
    n = f->n;
    t = f->t;
    A = f->terms;
    if( x<0 || x>=n ) return 0;
    m = termsize(n);
    q = x/64;
    r = x%64;
    mask = 1LU << r;
    //if( r==0 ) j = q; else j = q+1;  BUG?
    j = q;
    for( i=0; i<t; i++, A+=m )
        if( A[j]&mask ) return 1;
    return 0;
}


ULONG *degrees(int n, ULONG *X) {
    int i;
    ULONG *D;
    D = array64u(n);
    for( i=0; i<n; i++ ) 
        if( X[i]==0 ) D[i] = 0; else D[i] = 1L << degpoly(X[i]);
    return D;
}

void setbit(int x, ULONG *X) {
// set bit x of X to be 1
    int q,r; ULONG b;
    q = x/64;
    r = x%64;
    b = 1LU << r;
    X[q] = X[q] | b;
    return ;
}

void printhash(ULONG h,char x) {
    int i,first;
    if( h==0 ) { printf("0;\n"); return; }
    first = 1;
    for( i=0; h; i++ ) { 
        if( h&1 ) { 
            if( first ) first = 0; else printf("+");
            if( i==0 ) printf("1"); else printf("%c^%d",x,i);
        }
        h = h >> 1;
    }
    printf(";\n");
    return;
}


ULONG hashpoly(poly *f, ULONG *X, ULONG *B, ULONG minpoly) {
    int i,j,k,n,t,m;
    ULONG h,x,mask,*A;
    if( minpoly==0 ) { printf("hashpoly: minpoly=0\n"); exit(1); }
    mask = 1LU << degpoly(minpoly);
    n = f->n;
    t = f->t;
    A = f->terms;
    m = termsize(n);
    h = 0;
    for( i=0; i<t; i++,A+=m ) {
        x = hashprod(n,A,X,minpoly,B,mask);
        h = h ^ x; // h(z) = h(z) + x(z) mod 2
    }
    return h;
}

/**************************************************************************/
/*                                                                        */
/*                      Algebraic Library Here                            */
/*                                                                        */
/**************************************************************************/


int numxterms( poly *f, int x ) {
    // count the number of terms in f that have x
    int n,t,q,r,termsize,i,k; ULONG mask,*A;
    n = f->n;
    t = f->t;
    A = f->terms;
    if( x>=n ) { printf("evalto0: x is out of range\n"); exit(1); }
    q = n/64;
    r = n%64;
    if( r=0 ) termsize = q; else termsize = q+1;
    q = x/64;
    r = x%64;
    mask = 1;
    mask = mask << r;
    for( k=0, i=0; i<t; i++ ) { 
        if( A[q]&mask ) k++;
        A += termsize;
    }
    return k;
}


poly* xfactors( poly *f ) {
    // identify all trivial factors x_i that divide f(x1,...,xn)
    // return g(x1,...,xn) = product of all x_i|f, the monomial factor of f
    // In Maple: mul( x^ldegree(f,x), x in indets(f) );
    int n,t,m,i,j; ULONG *A, *X; 
    poly *g;
    n = f->n;
    t = f->t;
    A = f->terms;
    m = termsize(n);
    g = newpoly(n,1);
    X = g->terms;
    for( j=0; j<m; j++ ) X[j] = 0xFFFFFFFFFFFFFFFFL;
    for( i=0; i<t; i++, A+=m )
        for( j=0; j<m; j++ )
            X[j] = X[j] & A[j];
    return g;
}

poly* divxfactors( poly *f ) {
    // compute f/xfactors(f) to make f monomial primitive
    // e.g. f = x0*x1*x2+x0*x2*x3+x0*x3 --> x1*x2+x2*x3+x3
    int n,t,m,i,j; ULONG *A, *X; 
    poly *g,*h;
    n = f->n;
    t = f->t;
    A = f->terms;
    m = termsize(n);
    g = newpoly(n,1);
    X = g->terms;
    for( j=0; j<m; j++ ) X[j] = 0xFFFFFFFFFFFFFFFFL;
    for( i=0; i<t; i++, A+=m )
        for( j=0; j<m; j++ )
            X[j] = X[j] & A[j];
    for( j=0; j<m; j++ ) X[j] = ~X[j];
    h = copypoly(f);
    A = h->terms;
    for( i=0; i<t; i++, A+=m )
        for( j=0; j<m; j++ )
            A[j] = A[j] & X[j];
    freepoly(g);
    return h;
}

int isonepoly( poly *f ) {
    int i,n,m,t; ULONG *A;
    n = f->n;
    t = f->t;
    if( t != 1 ) return 0;
    m = termsize(n);
    A = f->terms;
    for( i=0; i<m; i++ ) if( A[i] ) return 0;
    return 1;
}

poly* evalxto0( poly *f, int x ) {
// eval(f,x=0) just delete terms with x in them
// Equivalent to coeff(f,x,0)
    int n,t,q,r,m,s,i,j,k;
    ULONG mask, *A, *B;
    poly *g;
    n = f->n;
    t = f->t;
    A = f->terms;
    q = n/64;
    r = n%64;
    if( r==0 ) m = q; else m = q+1;
    s = numxterms(f,x);
    B = array64u( m*(t-s) );
    g = (poly *) malloc(sizeof( poly ));
    g->n = n;
    g->t = t-s;
    g->terms = B;
    q = x/64;
    r = x%64;
    mask = 1;
    mask = mask << r;
    for( k=0,i=0; i<t; i++,A+=m ) {
        if( A[q]&mask ) continue;
        for( j=0; j<m; j++ ) B[j] = A[j];
        B+=m;
        k++;
    }
//printf("evalxto0: t=%d s=%d k=%d\n",t,s,k);
    return g;
}
    

poly* diffpoly( poly *f, int x ) {
// diff(f,x)
    int n,t,q,r,m,s,i,j,k;
    ULONG mask, *A, *B;
    poly *g;
    n = f->n;
    t = f->t;
    A = f->terms;
    q = n/64;
    r = n%64;
    if( r==0 ) m = q; else m = q+1;
    s = numxterms(f,x);
    B = array64u( m*s );
    g = (poly *) malloc(sizeof( poly ));
    g->n = n;
    g->t = s;
    g->terms = B;
    q = x/64;
    r = x%64;
    mask = 1;
    mask = mask << r;
    for( k=0,i=0; i<t; i++,A+=m ) {
        if( A[q]&mask ) {
            for( j=0; j<m; j++ ) B[j] = A[j];
            B[q] = B[q] ^ mask; // set degree to 0
            B+=m;
            k++;
        // else omit it
        }
    }
//printf("diffpoly: t=%d s=%d k=%d\n",t,s,k);
    return g;
}

poly* coeff( poly *f, int x ) { return diffpoly(f,x); }

int compare(int m, ULONG *X, ULONG *Y) {
// return 1 if X>Y, -1 if X<Y and 0 if X=Y
    int i;
    for( i=0; i<m; i++ ) {
       if( X[i]<Y[i] ) return -1;
       if( X[i]>Y[i] ) return 1;
    }
    return 0;
}

void sortswap(int m, ULONG *X, ULONG *Y) {
    int i; ULONG t;
    for( i=0; i<m; i++ ) { t = X[i]; X[i] = Y[i]; Y[i] = t; }
    return;
}

void sortmove( int m, ULONG *y, ULONG *x ) {
    int i;
    for( i=0; i<m; i++ ) y[i] = x[i];
    return;
}

void insertionsort(poly *f) {
// insertion sort on the terms in f
    int n,t,i,j,q,r,m;
    ULONG *A;
    n = f->n;
    t = f->t;
    A = f->terms;
    q = n/64;
    r = n%64;
    if( r==0 ) m = q; else m = q+1;
    for( i=1; i<t; i++ ) {
        for( j=i; j>0 && compare(m,A+j*m,A+(j-1)*m)==1; j-- )
            sortswap(m,A+j*m,A+(j-1)*m);
    }
    return;
}

int issorted(poly *f) {
    int n,t,i,j,q,r,m,c,flag;
    ULONG *A;
    n = f->n;
    t = f->t;
    A = f->terms;
    q = n/64;
    r = n%64;
    if( r==0 ) m = q; else m = q+1;
    //printf("issorted: ");
    flag = 1;
    for( i=1; i<t; i++ ) {
        c = compare(m,A+(i-1)*m,A+i*m);
        //if( c==-1 ) printf("<"); else if( c==1 ) printf(">"); else printf("=");
        if( c<0 ) flag = 0;
    }
    return flag;
}

int vec2partition64s(ULONG *D, int n, int b, ULONG *x, ULONG *y ) {
    int i,j,m;
    m = n/2;
    sortmove(b,x,D+m*b); sortmove(b,D+m*b,D+0);
    for( i=1,j=n-1; j>=i; ) {
       if( compare(b,D+i*b,x)>0 ) { sortmove(b,D+(i-1)*b,D+i*b); i++; }
       else { sortswap(b,D+i*b,D+j*b); j--; }
    }
    sortmove(b,D+(i-1)*b,x);
    return i-1;
}

void vec2sort64s(ULONG *A, int n, int b, ULONG *x, ULONG *y ) {
    int i;
    if( n<2 ) return;
    if( n>30 ) {  // use quicksort
        i = vec2partition64s(A,n,b,x,y);
        vec2sort64s(A,i,b,x,y);
        vec2sort64s(A+(i+1)*b,n-i-1,b,x,y);
    } else { // use insertion sort
        vec2sort64s(A,n-1,b,x,y);
        sortmove(b,x,A+(n-1)*b);
        for( i=n-1; i>0 && compare(b,x,A+(i-1)*b)==1; i-- ) sortmove(b,A+i*b,A+(i-1)*b);
        sortmove(b,A+i*b,x);
    }
    return;
}

void sortpoly(poly* f) 
{   // quicksort
    int n,m,t;
    ULONG *A,*x,*y;
    n = f->n;
    t = f->t;
    A = f->terms;
    m = termsize(n);
    x = array64u(m);
    y = array64u(m);
    vec2sort64s(A,t,m,x,y);
    free(x);
    free(y);
    return;
}

void mergesort2(int m, int t, ULONG *A, ULONG *B) {
    int s,i,j,k,r,c; ULONG *C;
    if( t<2 ) return;
    s = t/2;
    mergesort2(m,s,A,B);
    mergesort2(m,t-s,A+s,B+s);
    for( i=0; i<m*t; i++ ) B[i] = A[i];
    C = B+s*m;
    // merge B with C into A
    for( i=0,j=0; i<s && j<t-s; ) { 
        c = compare(m,B,C);
        if( c==1 ) { for( k=0; k<m; k++ ) A[k] = B[k]; A+=m; B+=m; i++; }
        else if( c==-1 ) { for( k=0; k<m; k++ ) A[k] = C[k]; A+=m; C+=m; j++; }
        else { for( k=0; k<m; k++ ) A[k] = B[k]; A+=m; B+=m; i++; 
               for( k=0; k<m; k++ ) A[k] = C[k]; A+=m; C+=m; j++; }
    }
    while(  i<s  ) { for( k=0; k<m; k++ ) A[k] = B[k]; A+=m; B+=m; i++; }
    while( j<t-s ) { for( k=0; k<m; k++ ) A[k] = C[k]; A+=m; C+=m; j++; }
    return;
}

void mergesort( poly *f ) {
    int n,t,m; ULONG *A,*B;
    n = f->n;
    t = f->t;
    if( t<50 ) return sortpoly(f);
    m = termsize(n);
    A = f->terms;
    B = array64u(t*m); // working storage
    mergesort2(m,t,A,B);
    free(B);
    return;
}


poly* project( poly *f, ULONG * X ) {
// Input f in n variables, X a bit vector of length n
// Output Y the projection of f onto the variables in X: Y[i]==1 iff haspoly(f,X[i])
    int n,t,m,i,j,k;
    ULONG *A, *B;
    poly *g,*h;
    n = f->n;
    t = f->t;
    A = f->terms;
    g = newpoly(n,t);
    B = g->terms;
    m = termsize(n);
    for( i=0; i<t; i++, A+=m, B+=m )
        for( j=0; j<m; j++ ) {
            B[j] = A[j] & X[j];
            //printf("A[%d]=%lu B[%d]=%lu\n",j,A[j],j,B[j]);
        }
    //sortpoly(g);
    mergesort(g);
    //if( issorted(g) ) printf("is sorted\n"); else printf("not sorted\n");
    A = g->terms;
    for( i=0,k=0; i<t; i=j,k++ ) {
        for( j=0; j<m; j++ ) A[k*m+j] = A[i*m+j];
        for( j=i+1; j<t && compare(m,A+m*i,A+m*j)==0; j++ ); 
    }
    if( k==t ) return g;
    h = newpoly(n,k);
    B = h->terms;
    //for( i=0,k=0; i<t; i=j,k++ ) {
    //   for( j=0; j<m; j++ ) B[k*m+j] = A[i*m+j];
    //   for( j=i+1; j<t && compare(m,A+m*i,A+m*j)==0; j++ );
    //}
    for( i=0; i<k; i++, A+=m, B+=m )
        for( j=0; j<m; j++ ) B[j] = A[j];
    freepoly(g);
    return h;
}

int equalpoly( poly *f, poly *g ) {
    int n,s,t,i,j,m; ULONG *A,*B;
    n = f->n;
    if( g->n != n ) { printf("addpoly: f and g must have same variables\n"); exit(1); }
    s = f->t;
    t = g->t;
    if( s!=t ) return 0;
    A = f->terms;
    B = g->terms;
    m = termsize(n);
    for( i=0; i<t; i++,A+=m,B+=m ) 
        for( j=0; j<m; j++ ) if( A[j]!=B[j] ) return 0;
    return 1;
}

poly* addpoly( poly *f, poly *g ) {
    int n,s,t,q,r,m,i,j,k,l,c;
    ULONG *A,*B,*C;
    poly *h;
    n = f->n;
    if( g->n != n ) { printf("addpoly: f and g must have same variables\n"); exit(1); }
    s = f->t;
    t = g->t;
    A = f->terms;
    B = g->terms;
    q = n/64;
    r = n%64;
    if( r==0 ) m = q; else m = q+1;
    h = newpoly(n,s+t);
    C = h->terms;
    i=0,j=0,k=0;
    while( i<s && j<t ) {
        c = compare(m,A+i*m,B+j*m);
        if( c==1 ) { for( l=0; l<m; l++ ) C[k*m+l] = A[i*m+l]; i++; k++; }
        if( c==0 ) { i++; j++; }
        if( c==-1 ) { for( l=0; l<m; l++ ) C[k*m+l] = B[j*m+l]; j++; k++; }
    }
    while( i<s ) { for( l=0; l<m; l++ ) C[k*m+l] = A[i*m+l]; i++; k++; }
    while( j<t ) { for( l=0; l<m; l++ ) C[k*m+l] = B[j*m+l]; j++; k++; }
    return h;
}

poly* mulpoly( poly *f, poly *g ) {
    // f(x0,x1,x2) x g(x3,x4,x5,x6) over F2
    int n,s,t,q,r,m,i,j,k,l;
    ULONG *A,*B,*C;
    poly *h;
    n = f->n;
    if( g->n != n ) { printf("mulpoly: f and g must have same variables\n"); exit(1); }
    s = f->t;
    t = g->t;
    A = f->terms;
    B = g->terms;
    q = n/64;
    r = n%64;
    if( r==0 ) m = q; else m = q+1;
    h = newpoly(n,s*t);
    C = h->terms;
    k = 0;
    for( j=0; j<t; j++ )
        for( i=0; i<s; i++,k++ )
            for( l=0; l<m; l++ ) C[k*m+l] = A[i*m+l] | B[j*m+l];
    sortpoly(h);
    return h;
}

poly* enlargepoly( poly *f ) {
    int n,t,i,j,k,m,s;
    ULONG x,y,z,*A,*B,mask;
    poly *g;
    n = f->n;
    t = f->t;
    A = f->terms;
    g = newpoly(2*n,t);
    m = termsize(n);
    B = g->terms;
    s = termsize(2*n);
    for( i=0; i<t; i++,A+=m,B+=s )
        for( j=0; j<m; j++ ) {
            x = A[j];
            for( y=0,k=0,mask=1; k<32; k++,mask<<=1 )
                y = y | ((x&mask) << k);
            B[2*j] = y;
            for( z=0,k=32; k<64; k++,mask<<=1 )
                z = z | ((x&mask) << k);
            //printf("enlarge: x=%lu y=%lu z=%lu\n",x,y,z);
            if( z ) B[2*j+1] = z;
        }
    B = g->terms;
    return g;
}

poly* mulpoly2( poly *f, poly *g ) {
    // f(x0,x1,x2) x g(x3,x4,x5,x6) allowing two bits for monomials 
    // the only difference is we add exponent vectors not xor them
    int n,s,t,q,r,m,i,j,k,l;
    ULONG *A,*B,*C;
    poly *h;
    n = f->n;
    if( g->n != n ) { printf("mulpoly: f and g must have same variables\n"); exit(1); }
    f = enlargepoly(f);
    g = enlargepoly(g);
    n = f->n;
    s = f->t;
    t = g->t;
    A = f->terms;
    B = g->terms;
    q = n/64;
    r = n%64;
    if( r==0 ) m = q; else m = q+1;
    h = newpoly(n,s*t);
    C = h->terms;
    k = 0;
    for( j=0; j<t; j++ )
        for( i=0; i<s; i++,k++ )
            for( l=0; l<m; l++ ) C[k*m+l] = A[i*m+l] + B[j*m+l];
    sortpoly(h);
    freepoly(f);
    freepoly(g);
    return h;
}

listpoly* facpoly( poly *f ) {
    int n,t,i,j,m,x;
    ULONG minpoly,*Z,z,d;
    ULONG a1,a2,b1,b2,c1,c2,d1,d2,h1,h2,*X1,*X2,*Y1,*Y2;
    poly *g,*h,*A,*B,*C,*D;
    listpoly *F,*G;

clock_t st,ht,et,dt,pt;

    n = f->n;
    t = f->t;
    if( DEBUGfacpoly ) if( t<20 ) { printf("facpoly: f:="); printpoly(f,'x'); }
    if( t<=1 ) { printpoly(f,'x'); F = conspoly(f,NULL); return F; }

st = clock();

    F = NULL; // list of factors

    g = xfactors(f); // the monomial content of f

    if( !isonepoly(g) ) {
        printf("facpoly: moncont(f)="); printpoly(g,'x'); 
        F = conspoly(g,F);
        f = divxfactors(f);
    }
//pt = clock()-st;
pt = 0;

    if( t<=3 ) { printpoly(f,'x'); F = conspoly(f,F); return F; }

ht = 0;
dt = 0;
counter = 0;

    // set up hashing into GF(2^63)
    minpoly = 16422940932436772585LU;
    d = degpoly(minpoly);
    X1 = randbits(d,n);
    Y1 = degrees(n,X1);
    X2 = randbits(d,n);
    Y2 = degrees(n,X2);
    m = termsize(n);

    // pick the first x such that deg(f,x)>0
    x = 0; while( x<n && !haspoly(f,x) ) x++;
    printf("facpoly: x=%d\n",x+1);

et = clock();
    A = evalxto0(f,x); //printf("A := "); printpoly(A,'x'); 
    B = diffpoly(f,x); //printf("B := "); printpoly(B,'x'); 
dt += clock()-et;

et = clock();
    a1 = hashpoly(A,X1,Y1,minpoly);
    a2 = hashpoly(A,X2,Y2,minpoly);
    b1 = hashpoly(B,X1,Y1,minpoly);
    b2 = hashpoly(B,X2,Y2,minpoly);
    //printf("a1=%lu a2=%lu b1=%lu, b2=%lu\n",a1,a2,b1,b2);
ht += clock()-et;

    Z = array64u(m);
    for( i=0,j=-1; i<n; i++ ) {
       if( i%64==0 ) { z = 0; j++; }
       if( i==x ) { h1 = 1; h2 = 1; }
       else {
et = clock();
          C = diffpoly(A,i); //printf("C := "); printpoly(C,'x');
          D = diffpoly(B,i); //printf("D := "); printpoly(D,'x');
dt += clock()-et;
et = clock();
          c1 = hashpoly(C,X1,Y1,minpoly);
          c2 = hashpoly(C,X2,Y2,minpoly);
          d1 = hashpoly(D,X1,Y1,minpoly);
          d2 = hashpoly(D,X2,Y2,minpoly);
          //printf("c1=%lu c2=%lu d1=%lu, d2=%lu\n",c1,c2,d1,d2);
          h1 = mulmod(a1,d1,minpoly) ^ mulmod(b1,c1,minpoly);
          h2 = mulmod(a2,d2,minpoly) ^ mulmod(b2,c2,minpoly);
ht += clock()-et;
          //printf("x%d: h1=%lu  h2=%lu  #C=%d  #D=%d\n",i,h1,h2,C->t,D->t);
          freepoly(C);
          freepoly(D);
       }
       if( h1 || h2 ) z = z | (1L << (i%64) );
       Z[j] = z;
    }

et = clock();
    g = project(f,Z); // g is irreducible
    for( j=0; j<m; j++ ) Z[j] = ~Z[j];
    h = project(f,Z);
pt += clock()-et;

if( g->t == 1 && h->t > 1 ) {
    printf("h = "); printpoly(h,'y');
    exit(1);
}

et = clock()-st;

    if( f->t > 500 ) { double etf,htf,cnt,dtf,ptf;
        etf = et/1000000.0; // convert to seconds
        htf = ht/1000000.0;
        htf = 100*htf/etf; // %age of time in GF(2^63)
        dtf = dt/1000000.0;
        dtf = 100*dtf/etf;
        ptf = pt/1000000.0;
        ptf = 100*ptf/etf;
        cnt = counter/1000000.0;
        printf("facpoly: time=%9.4f (%5.2f\%) (%5.2f\%) (%5.2f\%) #muls=%6.2f million\n",etf,htf,dtf,ptf,cnt);
    }

    t = g->t;
    if( DEBUGfacpoly )
        if( t<10 ) { printf("facpoly: g:="); printpoly(g,'x'); } else printf("facpoly: #g=%d\n",t);
    if( !isonepoly(g) ) F = conspoly(g,F);
    t = h->t;
    if( DEBUGfacpoly )
        if( t<10 ) { printf("facpoly: h:="); printpoly(h,'x'); } else printf("facpoly: #h=%d\n",t);
    if( !isonepoly(h) ) {
        G = facpoly(h);
        F = joinpoly(F,G);
    }

    return F;
}


listpoly* facpoly2( poly *f ) {
// Use the FD explicit algorithm
    int n,t,i,j,m,x,hash;
    ULONG z,*Z;
    poly *g,*h,*A,*B,*C,*D,*AD,*BC;
    clock_t st,ht,et;
    listpoly *F;

    n = f->n;
    t = f->t;
    if( t<10 ) { printf("facpoly2: f:="); printpoly(f,'x'); }
    if( t<=1 ) { printpoly(f,'x'); F = conspoly(f,NULL); return F; }

    F = NULL; // list of factors

    g = xfactors(f); // the monomial content of f

    if( !isonepoly(g) ) {
        printf("facpoly2: moncont(f)=");
        printpoly(g,'x'); 
        F = conspoly(g,F);
        f = divxfactors(f);
    }

    if( t<=3 ) { printpoly(f,'x'); F = conspoly(f,F); return F; }

st = clock();
ht = 0;

    // pick the first x such that deg(f,x)>0
    x = 0; while( x<n && !haspoly(f,x) ) x++;
    printf("facpoly2: x=%d\n",x);

    A = evalxto0(f,x); 
    B = diffpoly(f,x); 
//printf("facpoly2: A="); printpoly(A,'x');
//printf("facpoly2: B="); printpoly(B,'x');

    m = termsize(n);
    Z = array64u(m);
    for( i=0,j=-1; i<n; i++ ) {
       if( i%64==0 ) { z = 0; j++; }
       if( i==x ) hash = 1;
       else {
          C = diffpoly(A,i);
          D = diffpoly(B,i); 
//printf("C = "); printpoly(C,'x');
//printf("D = "); printpoly(D,'x');
et = clock();
          AD = mulpoly2(A,D);
          BC = mulpoly2(B,C);
//printf("AD = "); printpoly2(AD,'x');
//printf("BC = "); printpoly2(BC,'x');
          hash = !equalpoly(AD,BC);
          printf("#AD=%d  #BD=%d  AD=BC? %d\n",AD->t,BC->t,hash);
ht += clock()-et;
          freepoly(C);
          freepoly(D);
          freepoly(AD);
          freepoly(BC);
       }
       if( hash!=0 ) z = z | (1L << (i%64) );
       Z[j] = z;
    }

    g = project(f,Z); // g is irreducible
    for( j=0; j<m; j++ ) Z[j] = ~Z[j];
    h = project(f,Z);

et = clock()-st;

    t = g->t;
    if( t<10 ) { printf("facpoly2: g:="); printpoly(g,'x'); }
    else printf("facpoly2: #g=%d\n",t);
    if( !isonepoly(g) ) F = conspoly(g,F);
    t = h->t;
    if( t<10 ) { printf("facpoly2: h:="); printpoly(h,'x'); }
    else printf("facpoly2: #h=%d\n",t);
    if( !isonepoly(h) ) F = conspoly(h,F);

    if( f->t > 100 ) { double etf,htf,cnt,ptf;
        etf = et/1000000.0; // convert to seconds
        htf = ht/1000000.0;
        htf = 100*htf/etf; // %age of time in GF(2^63)
        printf("facpoly2: time=%9.4f (%5.2f\%)\n",etf,htf);
    }
    return F;
}

/*

int main() {
int i,j,k,n,s,t;
ULONG a,b,c,m,M,B,*A,*X,*D;
poly *f,*g,*h;
listpoly *F;

    a = 13;    // 1101;
    b = 5;     // 101;
    B = 4;     // 100;
    m = 23;    // 10111; 
    M = 16;    // 10000;
    extern ULONG seed,mult;
    seed = 1;
    mult = 6364136223846793003ll;
    counter = 0;
    c = MULMOD(a,b,m,B,M);
    printf("a=%u b=%u c=%u\n",a,b,c);
    b = 7;     // 111
    c = MULMOD(a,b,m,B,M);
    printf("a=%u b=%u c=%u\n",a,b,c);

    f = randpoly(10,9);
    printf("f = "); printpoly(f,'x');
    for( i=0; i<3; i++ ) printf("#x%d = %d\n",i,numxterms(f,i));
    g = evalxto0(f,2);
    printf("eval(f,x2=0) = "); printpoly(g,'x');
    g = diffpoly(f,2);
    printf("diff(f,x2) = "); printpoly(g,'x');
    M = 6;
    g = project(f,&M);
    printf("project(f,x1x2) = "); printpoly(g,'x');
    printf("f = "); printpoly(f,'x');
    printf("sorted(f): %d\n",issorted(f));
    //sortpoly(f);
    //printf("sorted(f): %d\n",issorted(f));
    printf("f = "); printpoly(f,'x');
    f = newpoly(4,4);
    A = f->terms;
    A[0] = 13;
    A[1] = 9;
    A[2] = 15;
    A[3] = 10;
    printf("f = "); printpoly(f,'x');
    g = xfactors(f);
    printf("g = "); printpoly(g,'x');
    h = divxfactors(f);
    printf("h = "); printpoly(h,'x');
    printf("factoring f = "); printpoly(f,'x');

facpoly(f);
printf("factor2 starting\n");
facpoly2(f);


    f = newpoly(4,4);
    A = f->terms;
    A[0] = 5;
    A[1] = 9;
    A[2] = 6;
    A[3] = 10;
    // f = x0*x1 + x0*x3 + x1*x2 + x2*x4;
    printf("f = "); printpoly(f,'x');
    X = indetspoly(f);
    printf("indets(f)=%lu\n",X[0]);
    // f2 = x0*x1 + x0*x1^2 + x0^2*x1 + x0^2*x2^2
    printf("f2 = "); printpoly2(f,'x');

facpoly(f);
printf("\nfactor2 starting\n");
facpoly2(f);


    printf("f = "); printpoly(f,'x');
    g = enlargepoly(f);
    printf("enlarge(f) = "); printpoly(g,'x');

    M = 12;
    g = project(f,&M); // project onto x2,x3
    printf("g = proj(f,x2x3) = "); printpoly(g,'x');
    M = 3;
    h = project(f,&M); // project onto x0,x1
    printf("h = proj(f,x0x1) = "); printpoly(h,'x');
    f = mulpoly( h, g );

facpoly(f);
printf("\nfactor2 starting\n");
facpoly2(f);

    printf("g = "); printpoly(g,'x');
    printf("h = "); printpoly(h,'x');

    printf("g^2 = "); printpoly2(mulpoly2(g,g),'x');
    printf("f^2 = "); printpoly2(mulpoly2(f,f),'x');

    printf("f = h x g = "); printpoly(f,'x');
    h = randpoly1(6,3,5);
    printf("randpoly1: #=%d ",h->t); printpoly(h,'y');
    h = randpoly2(8,4,6);
    printf("randpoly2: #=%d ",h->t); printpoly(h,'y');
    h = randpoly2(10,5,6);
    printf("randpoly2: #=%d ",h->t); printpoly(h,'y');

    h = randpoly1(180,10,2);
    printf("randpoly1: h= "); printpoly(h,'x');
    shiftpoly(h,10);
    printf("shiftpoly: h= "); printpoly(h,'x');
    shiftpoly(h,54);
    printf("shiftpoly: h= "); printpoly(h,'x');
    shiftpoly(h,10);
    printf("shiftpoly: h= "); printpoly(h,'x');
    shiftpoly(h,60);
    printf("shiftpoly: h= "); printpoly(h,'x');
    h = randpoly1(180,10,2);
    printf("randpoly1: h= "); printpoly(h,'x');
    shiftpoly(h,100);
    printf("randpoly1: h= "); printpoly(h,'x');

{
int k;
poly *A,*B,*C,*D,*AD,*BC;
ULONG a,b,c,d,x,minpoly,*X,*Y,*Z;
    A = evalxto0(f,0);
    B = diffpoly(f,0);
    printf("f="); printpoly(f,'x');
    printf("A="); printpoly(A,'x');
    printf("B="); printpoly(B,'x');
    printf("\n");
    printf("f := "); printpoly(f,'x');
    minpoly = 64+32+1; // z^6+z^5+1
    d = degpoly(minpoly);
    printf("m := ",minpoly); printhash(minpoly,'z');
    n = f->n;
    X = randbits(d,n);
    Y = degrees(n,X);
    for( i=0; i<n; i++ ) { printf("x%d := ",i); printhash(X[i],'z'); }
    x = hashpoly(f, X, Y, minpoly);
    printf("#x = %lu = ",x); printhash(x,'z');
    printf("\n");
    printf("A := "); printpoly(A,'x');
    a = hashpoly(A,X,Y,minpoly); printf("#a = %lu = ",a); printhash(a,'z');
    printf("B := "); printpoly(B,'x');
    b = hashpoly(B,X,Y,minpoly); printf("#b = %lu = ",b); printhash(b,'z');

    C = diffpoly(A,1);
    D = diffpoly(B,1);
    printf("C := "); printpoly(C,'x');
    c = hashpoly(C,X,Y,minpoly); printf("#c = %lu = ",c); printhash(c,'z');
    printf("D := "); printpoly(D,'x');
    d = hashpoly(D,X,Y,minpoly); printf("#d = %lu = ",d); printhash(d,'z');
    printf("a*d=%lu\n",mulmod(a,d,minpoly));
    printf("b*c=%lu\n",mulmod(b,c,minpoly));
    x = mulmod(a,d,minpoly) ^ mulmod(b,c,minpoly);
    printf("#AD-BC=%lu\n",x); printf("x := "); printhash(x,'z');
    printf("\n");

    printf("\n");
    AD = mulpoly2(A,D);
    BC = mulpoly2(B,C);
    printf("AD = "); printpoly2(AD,'x');
    printf("BC = "); printpoly2(BC,'x');
    if( equalpoly(AD,BC) ) printf("AD=BC\n"); else printf("AD<>BC\n");
    printf("\n");

    C = diffpoly(A,2);
    D = diffpoly(B,2);
    printf("C := "); printpoly(C,'x');
    printf("D := "); printpoly(D,'x');
    c = hashpoly(C,X,Y,minpoly); printf("#c = %lu = ",c); printhash(c,'z');
    d = hashpoly(D,X,Y,minpoly); printf("#d = %lu = ",d); printhash(d,'z');
    x = mulmod(a,d,minpoly) ^ mulmod(b,c,minpoly);
    printf("#AD-BC=%lu\n",x); printf("H := "); printhash(x,'z');
    printf("\n");

    AD = mulpoly2(A,D);
    BC = mulpoly2(B,C);
    printf("AD = "); printpoly2(AD,'x');
    printf("BC = "); printpoly2(BC,'x');
    if( equalpoly(AD,BC) ) printf("AD=BC\n"); else printf("AD<>BC\n");
    printf("\n");

    C = diffpoly(A,3);
    D = diffpoly(B,3);
    printf("C := "); printpoly(C,'x');
    printf("D := "); printpoly(D,'x');
    c = hashpoly(C,X,Y,minpoly); printf("#c = %lu = ",c); printhash(c,'z');
    d = hashpoly(D,X,Y,minpoly); printf("#d = %lu = ",d); printhash(d,'z');
    x = mulmod(a,d,minpoly) ^ mulmod(b,c,minpoly);
    printf("#AD-BC=%lu\n",x); printf("H := "); printhash(x,'z');

    n = 5;
    f = newpoly(n,4);
    Z = f->terms;
    Z[0] = 3;
    Z[1] = 5;
    Z[2] = 6;
    Z[3] = 7;
    printf("f = "); printpoly(f,'x');
    sortpoly(f);
    printf("f = "); printpoly(f,'x');
    g = newpoly(n,3);
    Z = g->terms;
    Z[0] = 8;
    Z[1] = 16;
    Z[2] = 0;
    printf("g = "); printpoly(g,'x');
    sortpoly(g);
    printf("g = "); printpoly(g,'x');
    h = mulpoly(f,g);
    printf("h = "); printpoly(h,'x');
    A = evalxto0(h,0);
    B = diffpoly(h,0);
    printf("A="); printpoly(A,'x');
    printf("B="); printpoly(B,'x');
    printf("\n");
    for( i=1; i<n; i++ ) {
        C = diffpoly(A,i);
        D = diffpoly(B,i);
        a = hashpoly(A,X,Y,minpoly);
        b = hashpoly(B,X,Y,minpoly);
        c = hashpoly(C,X,Y,minpoly);
        d = hashpoly(D,X,Y,minpoly);
        printf("a=%lu b=%lu c=%lu d=%lu\n",a,b,c,d);
        x = mulmod(a,d,minpoly) ^ mulmod(b,c,minpoly);
        printf("x%i: h=%lu\n",i,x);
    }
    printf("\n");

    minpoly = 16422940932436772585LU;
    printf("m := ",minpoly); printhash(minpoly,'z');
    d = degpoly(minpoly);
    X = randbits(d,n);
    Y = degrees(n,X);
    for( i=1; i<n; i++ ) {
        C = diffpoly(A,i);
        D = diffpoly(B,i);
        a = hashpoly(A,X,Y,minpoly);
        b = hashpoly(B,X,Y,minpoly);
        c = hashpoly(C,X,Y,minpoly);
        d = hashpoly(D,X,Y,minpoly);
        x = mulmod(a,d,minpoly) ^ mulmod(b,c,minpoly);
        printf("x%i: h=%lu\n",i,x);
    }
    n = 10;
    t = 10; 
    d = degpoly(minpoly);
    X = randbits(d,n);
    Y = degrees(n,X);
    f = randpoly1( n, n/2, t );
    g = randpoly2( n, n/2, t );
    printf("f:="); printpoly(f,'x');
    printf("g:="); printpoly(g,'x');
    h = mulpoly(f,g);

printf("\n");
F = facpoly(h);
printf("Factors found\n");
while( F!=NULL ) { printf("  #f=%d",F->f->t); F = F->next; }
printf("\n");

    //printf("h:="); printpoly(h,'x');
    printf("#f=%d #g=%d #h=%d\n",f->n,g->n,h->n);
    A = evalxto0(h,0);
    B = diffpoly(h,0);
    printf("A:="); printpoly(A,'x');
    printf("B:="); printpoly(B,'x');
    printf("\n");
    for( i=1; i<n; i++ ) {
        C = diffpoly(A,i);
        D = diffpoly(B,i);
        a = hashpoly(A,X,Y,minpoly);
        b = hashpoly(B,X,Y,minpoly);
        c = hashpoly(C,X,Y,minpoly);
        d = hashpoly(D,X,Y,minpoly);
        x = mulmod(a,d,minpoly) ^ mulmod(b,c,minpoly);
        printf("x%i: h=%lu\n",i,x);
    }

  int nn[3] = {20,50,80};
  //int ss[4] = {10,25,50,100};
  //int tt[4] = {100,40,20,10};
  int ss[4] = {10,25,50,100};
  int tt[4] = {1000,400,200,100};
  //int ss[4] = {50,100,200,316};
  //int tt[4] = {2000,1000,500,316};
  //int ss[4] = {100,200,500,1000};
  //int tt[4] = {10000,5000,2000,1000};
  for( i=0; i<1; i++ ) {
    n = 100;
    j = nn[i];
    for( k=0; k<4; k++ ) {
      s = ss[k]; t = tt[k];
      d = degpoly(minpoly);
      X = randbits(d,n);
      Y = degrees(n,X);
      f = randpoly1( n, j, s );
      g = randpoly2( n, j, t );
      //printf("f:="); printpoly(f,'x');
      //printf("g:="); printpoly(g,'x');
      h = mulpoly(f,g);
      //printf("h:="); printpoly(h,'x');
      printf("#f=%d #x=%d #g=%d #h=%d\n",f->t,j,g->t,h->t);
      //F = facpoly2(h);
      F = facpoly(h);
      printf("Factors found\n");
      while( F!=NULL ) { printf("  #f=%d",F->f->t); F = F->next; }
      printf("\n");
    }
    printf("\n\n");
  }
}
    return 1;
}

*/