// IPMUL6.c
// Add IGCDCOF
// Copyright, Michael Monagan, April 2022, May 2022, February 2023
// Compile with:  gcc -O3 polyalg8.so inv1.so -shared -o IPMUL.so -fPIC IPMUL7.c

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

#define debug 0
#define LONG long long int

// In the following routines R is the quotient ring Zp[x1,x2,...,xn]/<m1,m2,...,mn>
// where p is a prime < 2^63 stored as a LONG integer.
// m1 is in R2[x1] where R2=Zp[x2,...,xn]/<m2,...,mn>,
// m2 is in R3[x2] where R3=Zp[x3,...,xn]/<m3,...,mn>, ... and 
// mn is in Zp[xn] so that T = {m1,m2,...,mn} forms a triangular set.
// The quotient ring R is isomorphic to the vector space Zp^d of dimension d.
// Here (and in the codes below) d = product( degree(m_i,x_i), 1<=i<=n ).
//
// The "minimal polynomials" m1,m2,...,mn are monic in xi but not necessarily 
// irreducible so that R is really a finite ring, not necessarily a field.
// They are encoded in the array M below and D is an array of their degrees.
// The number of minimal polynomials N and their degrees in D define the 
// space needed for them and for elements of R.
//
// Most routines take N,D,M,p as inputs.
//
// The polynomials are encoded densely in an array with the degree stored.
// A zero polynomial has degree -1 by convention.
// For example, for N=1, m=y^3-2, D=[3], a = 3*y*x^2+4*y^2+5, we encode a as
// [deg(a,x),coeff(a,x,0),coeff(a,x,1),...] = [2,(2,5,0,4),(-1,0,0,0),(1,0,3,0)]
// Since D=[3] each coefficient is in R=Zp[y]/<m> and has at most degree 2.
//
// The exported (public) routines are the following.
// Here a,b,c,s,gamma are in R and A,B,C are in R[x]

void IPADD( LONG *A, LONG *B, int N, int *D, LONG p ); // A += B
void IPSUB( LONG *A, LONG *B, int N, int *D, LONG p ); // A -= B
void IPREM( LONG *A, LONG *B, int N, int *D, LONG *M, LONG *T, LONG p ); // A = A rem B
void IPMUL( LONG *A, LONG *B, LONG *C, int N, int *D, LONG *M, LONG *T, LONG p ); // C = A B
void IPELTMUL( LONG *a, LONG *b, LONG *c, int N, int *D, LONG *M, LONG *T, LONG p ); // c = a b
void IPSCAMUL( LONG *s, LONG *A, int N, int *D, LONG *M, LONG *T, LONG p ); // A = s A
int IPINV( LONG *A, int N, int *D, LONG *M, LONG *I, LONG *W, LONG p ); // I = A^(-1)
int IPGCD( LONG *A, LONG *B, int N, int *D, LONG *M, LONG *Z, LONG *W, LONG p ); // A = gcd(A,B)
int IPGCDCOF( LONG *A, LONG *B, LONG *CA, LONG *CB, int N, int *D, LONG *M, LONG *Z, LONG *W, LONG p );
int PHIGCD( LONG *A, LONG *B, int N, int *D, LONG *M, LONG *gamma, LONG *Z, LONG *W, LONG p );
int PHIGCDCOF( LONG *A, LONG *B, LONG *CA, LONG *CB, int N, int *D, LONG *M, LONG *gamma, LONG *Z, LONG *W, LONG p );
int PHICOF2( LONG *A, LONG *B, LONG *G, int N, int *D, LONG *M, LONG *gamma, LONG *Z, LONG *W, LONG p );
int MONIC( LONG *A, int N, int *D, LONG *M, LONG *Z, LONG *W, LONG p ); // A = A / LC(A)

// These are conversions from REDCEN to DENSE
void IPremovedegreesPOLY( LONG *A, LONG *B, int N, int *D ); // B may equal A
void IPremovedegreesRING( LONG *A, LONG *B, int N, int *D ); // B may equal A
void IPinsertdegreesPOLY( LONG *A, int d, LONG *B, int N, int *D ); // B may equal A
void IPinsertdegreesRING( LONG *A, LONG *B, int N, int *D ); // B may equal A

// These routines are defined in polyalg7.c

LONG mod64s(LONG a, LONG p);  // a mod p
LONG mul64s(LONG a, LONG b, LONG p);  // a*b mod p
LONG modinv64s(LONG a, LONG p); // a^(-1) mod p
LONG * array(LONG n);
int * arrayint(LONG n);
void vecprint64s( LONG *A, int n );
void veccopy64s( LONG *u, int n, LONG *v );
void polprint64s( LONG *A, int d, LONG p );
void polcopy64s( LONG *A, int d, LONG *B );
int polmul64s( LONG * A, LONG * B, LONG * C, int da, int db, LONG p);
int poldiv64s( LONG * A, LONG * B, int da, int db, LONG p );
int polgcd64s( LONG * A, LONG * B, int da, int db, LONG p );
int poladdto64s( LONG *A, LONG *B, int da, int db, LONG p );
int polsubfrom64s( LONG *A, LONG *B, int da, int db, LONG p );
// Solve S A + T B = G in Zp[x]
void polgcdext64s( LONG *A, LONG *B, int da, int db,
                   LONG *G, LONG *S, LONG *T, int *dG, int *dS, int *dT,
                   LONG *W, LONG p );

// These are utilities for computing space requirements
int SpaceElem( int N, int *D ); // to store an element of R
int SpacePoly( int d, int N, int *D ); // to store f in R[x] of degree d
int SpaceMinpoly( int N, int *D ); // to store the main minimal polynomial in R
int TempSpace( int N, int *D ); // for multiplication in R
int TempSpaceIPMUL( int N, int *D ); // for multipication and division in R[x]
int TempSpaceIPINV( int N, int *D ); // for an inverse in R
int TempSpaceIPGCD( int N, int *D ); // for gcd in R[x] -- need inverse in R


/************************************************************************************/
/*******************************      Utilites      *********************************/
/************************************************************************************/

int MIN(int a, int b) { if(a<b) return a; else return b; } // inline
int MAX(int a, int b) { if(a>b) return a; else return b; } // inline
long NEG(LONG a,LONG p) { if( a==0 ) return 0; else return p-a; } // inline


void printpol( LONG *A, int N, int *D, char *X ) {
    int i,da,sp;
    da = A[0];
    if( da==-1 ) { printf("0"); return; }
    char x = X[0];
    A++;
    if( N==0 ) {
         for( i=0; i<=da; i++ ) {
             printf("%lld",A[i]);
             if( i>0 ) printf("*%c",x);
             if( i>1 ) printf("^%d",i);
             if( i<da ) printf("+");
         }
         return;
    }
    //sp = 1; for( i=N-1; i>=0; i++ ) sp = sp*(D[i]+1);
    sp = SpaceElem(N,D);
    for( i=0; i<=da; i++ ) {
        printf("(");  
        printpol(A+i*sp,N-1,D+1,X+1);
        printf(")");
        if(i>0) printf("*%c",x);
        if(i>1) printf("^%d",i);
        if(i<da) { if( N>1 ) printf(" + "); else printf("+"); }
    }
}


int SpaceElem( int N, int *D ) {
    // Compute space of one field element.  I had lots of errors in this.
    // It's not symmetric in counting up to N-1 or down from N-1
    // The degree of the first extension is D[N-1] not D[0] note
    int i,S;
    for( S=1,i=N-1; i>=0; i-- ) S = D[i]*S + 1;
    return S;
}

int SpacePoly( int d, int N, int *D ) {
    // Space for a polynomial of degree d in R[x]
    return (d+1)*SpaceElem(N,D)+1; // +1 for degree
}

int SpaceMinpoly( int N, int *D ) {
    return (D[0]+1)*SpaceElem(N-1,D+1)+1; // +1 for degree
}

int TempSpace( int N, int *D ) {
    // Compute temporary space for products of quotient ring elements
    // E.g R = Zp[x,y,z]/(x^3-yz,y^2-z-1,z^4-2) so that N=3 and D=[3,2,4]
    // A product in x has degree at most 4 in x, 1 in y and 3 in z so needs
    // (4+1)(2(3+1)+1))+1 = 46  where the +1's are for storing degrees
    // But we need to multiply this by 2 for two temporaries T1 and T2
    // and add the temporary space for N=2,D=[2,4] and for N=1,D=[4]
    int S,Sn,i,max_prod_deg;
    S = 0;
    Sn = 1;
    for( i=N-1; i>=0; i-- ) {
        max_prod_deg = 2*D[i]-1;
        S += (max_prod_deg*Sn+1)*2; // *2 for T1 and T2, 
        Sn = D[i]*Sn+1; // +1 for storing the degree
    }
    return S;
}

void COPYAtoB( LONG *A, LONG *B, int N, int *D ) {
    // CopyAtoB copies A into B where A is in R and B = array( SpaceElem(N,D) )
    int da,i,sp;
    if( N<1 ) return;
    B[0] = da = A[0];
    if( da==-1 ) return;
    A++; B++;
    if( N==1 ) { for( i=0; i<=da; i++ ) B[i] = A[i]; return; }
    sp = SpaceElem(N-1,D+1);
    for( i=0; i<=da; i++ ) COPYAtoB( A+i*sp, B+i*sp, N-1, D+1 );
    return;
}

void POLCOPYAtoB( LONG *A, LONG *B, int N, int *D ) {
    // POLCOPYAtoB copies A into B where A is in R[x]
    // B must be an array of size (da+1) SpaceElem(N,D) + 1
    int da,i,sp;
    B[0] = da = A[0];
    if( da==-1 ) return;
    A++; B++; // point to the coefficients
    if( N==0 ) { for( i=0; i<=da; i++ ) B[i] = A[i]; return; }
    sp = SpaceElem(N,D);
    for( i=0; i<=da; i++ ) COPYAtoB( A+i*sp, B+i*sp, N, D );
    return;
}

int MONIC( LONG *A, int N, int *D, LONG *M, LONG *Z, LONG *W, LONG p ) {
    // Make A in R[x] monic.  
    // Output the inverse of the leading coefficient in Z
    // A zero divisor is returned in Z.
    int da,spaceN,z,i,constant; LONG *LC;
    da = A[0]; if( da==-1 ) return 0;
    spaceN = SpaceElem(N,D);
    LC = A+1 + da*spaceN; // LC points to lcoeff(A,x)
    for( constant=1,i=0; i<N; i++ ) constant = constant && LC[i]==0;
    if( constant && LC[N]==1 ) {
        for( i=0; i<N; i++ ) Z[i] = 0;
        Z[N] = 1; // Z = 1 is the inverse
        return 0; // A is already monic
    }
    z = IPINV(LC,N,D,M,Z,W,p);
    if( z>0 ) return z;
    A[0] = da-1; // to avoid multiplying Z x LC = 1
    IPSCAMUL(Z,A,N,D,M,W,p);
    for( i=0; i<N; i++ ) LC[i]=0;
    LC[N] = 1; // A[da] = 1
    A[0] = da;
    return 0;
}

/***********************  EXPORTED  ROUTINES  *********************************/

void IPADD( LONG *A, LONG *B, int N, int *D, LONG p ) {
// A += B  where A,B are in R[x] - no temporary space is needed
    int da,db,i,sp,min;
    db = B[0];
    if( db==-1 ) return; // B = 0
    da = A[0];
    A++; B++;
    if( N==0 ) { 
        A[-1] = poladdto64s(A,B,da,db,p);  // A += B
        return;
    }
    sp = SpaceElem(N,D);
    min = MIN(da,db);
    for( i=0; i<=min; i++ )
        IPADD( A+i*sp, B+i*sp, N-1, D+1, p );
    for( ; i<=db; i++ )
        COPYAtoB( B+i*sp, A+i*sp, N, D );
    for( da=MAX(da,db); da>=0 && A[da*sp]==-1; da-- );
    A[-1] = da; // set the degree
    return;
}

void IPSUB( LONG *A, LONG *B, int N, int *D, LONG p ) {
// A -= B  where A,B are in R[x] - no temporary space is needed
    int da,db,i,sp,min;
    db = B[0];
    if( db==-1 ) return;
    da = A[0];
    A++; B++;
    if( N==0 ) {
        A[-1] = polsubfrom64s(A,B,da,db,p); // A -= B
        return;
    }
    sp = SpaceElem(N,D);
    min = MIN(da,db);
    for( i=0; i<=min; i++ )
        IPSUB( A+i*sp, B+i*sp, N-1, D+1, p );
    for( ; i<=db; i++ ) {
        A[i*sp] = -1;
        IPSUB( A+i*sp, B+i*sp, N-1, D+1, p );
    }
    for( da=MAX(da,db); da>=0 && A[da*sp]==-1; da-- );
    A[-1] = da; // set the degree
    return;
}

void IPCONMULrec( LONG c, LONG *A, LONG *B, int N, int *D, LONG p ) {
    int i,da,S,spaceN;
    da = A[0];
    B[0] = da;
    if( da<0 ) return;
    A++;
    B++;
    if( N==1 ) {
        if( c==1 ) for( i=0; i<=da; i++ ) B[i] = A[i];
        else for( i=0; i<=da; i++ ) B[i] = mul64s(c,A[i],p);
    } else {
        S = SpaceElem(N-1,D+1);
        for( i=0; i<=da; i++ ) IPCONMULrec(c,A+i*S,B+i*S,N-1,D+1,p);
    }
    return;
}

void IPCONMUL( LONG c, LONG *A, LONG *B, int N, int *D, LONG p ) {
    // B = c A where A in in R and c is a constant in Zp 
    // B may be the same array as A
    if( c==0 ) { B[0] = -1; return; }
    IPCONMULrec(c,A,B,N,D,p);
    return;
}

int ISCONSTANT( LONG *A, int N ) {
// 0 is not considered a constant note
    int i;
    for( i=0; i<N; i++ ) if( A[i]!=0 ) return 0;
    return 1;
}

void IPELTMUL( LONG *A, LONG *B, LONG *C, int N, int *D, LONG *M, LONG *T, LONG p ) {
    // C = A B where A, B are elements in R.
    // C = array( SpaceElem(N,T) ) and T = array( TempSpace(N,D) )
    int da,db,i,S,spaceM,spaceP;
    da = A[0];
    db = B[0];
    if( da==-1 || db==-1 ) { C[0] = -1; return; }
    if( N<=0 ) { printf("IPELEMMUL: N must be positive\n"); exit(1); }
    // E.g. for mutiplication of gamma x B where gamma = (c0 + c1 z1 + ... + ck zk)
    //if( ISCONSTANT(A,N) ) { IPCONMUL(A[N],B,C,N,D,p); return; }
    //if( ISCONSTANT(B,N) ) { IPCONMUL(B[N],A,C,N,D,p); return; }
    if( N==1 ) {
        T[0] = polmul64s(A+1,B+1,T+1,da,db,p);
        C[0] = poldiv64s(T+1,M+1,T[0],M[0],p);
        polcopy64s(T+1,C[0],C+1);
        return;
    }
    S = SpaceElem(N-1,D+1);
    spaceM = (D[0]+1)*S+1;
    spaceP = (2*D[0]-1)*S+1;
    IPMUL(A,B,T,N-1,D+1,M+spaceM,T+spaceP,p);
    IPREM(T,M,N-1,D+1,M+spaceM,T+spaceP,p);
    COPYAtoB(T,C,N,D);
    return;
}

void IPSCAMUL( LONG *s, LONG *A, int N, int *D, LONG *M, LONG *T, LONG p ) {
    // A = s A where s is a scalar in R and A is in R[x]
    int da,ds,i,S,spaceN,spaceM,spaceP;
    LONG *x;
    if( N==0 ) { printf("IPSCAMUL: should not happen\n"); exit(1); }
    ds = s[0]; if( ds==-1 ) { A[0] = -1; return; }
    da = A[0]; if( da==-1 ) return;
    for( S=1,i=N-1; i>0; i-- ) S = D[i]*S + 1;
    spaceM = (D[i]+1)*S+1;
    spaceP = (2*D[i]-1)*S+1;
    spaceN = D[i]*S+1;
    A++; // move to the coefficients!
    for( x=A,i=0; i<=da; i++,x+=spaceN ) {
        if( N==1 ) {
            T[0] = polmul64s(s+1,x+1,T+1,ds,x[0],p);
            x[0] = poldiv64s(T+1,M+1,T[0],M[0],p);
            polcopy64s(T+1,x[0],x+1);
        } else {
            IPMUL(s,x,T,N-1,D+1,M+spaceM,T+spaceP,p);
            IPREM(T,M,N-1,D+1,M+spaceM,T+spaceP,p);
            COPYAtoB(T,x,N,D);
        }
    }
    // s could be a zero divisor so s A[da] = 0 is possible
    while( da>=0 && A[da*spaceN]==-1 ) da--;
    A[-1] = da; // set degree
    return;
}

int TempSpaceIPMUL( int N, int *D ) {
    // space needed for all temporary products in IPMUL and IPREM
    int S,sp,i;
    S = 0;
    for( sp=1,i=N-1; i>=0; i-- ) {
        sp = D[i]*sp + 1;
        S += 2*sp; // two temporaries
    }
    return S;
}

int REDUCE( LONG *A, LONG z, LONG p ) {
// A = ax^2+bx+c  subs(x^2=-z,A)  ==>  A = bx+(c-az)
    int d; LONG t;
    A++;
    t = A[0]-mul64s(A[2],z,p);
    t += (t>>63) & p;
    A[0] = t;
    for( d=1; A[d]==0; d-- ); 
    return A[-1]=d;
}
 
void IPMUL( LONG *A, LONG *B, LONG *C, int N, int *D, LONG *M, LONG *T, LONG p ) {
    // Compute C = A B where A,B are in R[x] and C is space for the product:
    // C = array( ( deg(A)+deg(B)+1 ) SpaceElem(N,D) + 1 )
    // T = array( TempSpace(N,D) ) is temporary space for products in R
    int da,db,dc,dr,i,j,k,spaceN,spaceM,spaceP,bot,top; LONG z,*T1,*T2;
    da = A[0];
    db = B[0];
    if( da==-1 || db==-1 ) { C[0]=-1; return; }
    A++; B++; C++; // point to the coefficients
    if( N==0 ) { C[-1] = polmul64s(A,B,C,da,db,p); return; }
    spaceN = 1;
    for( i=N-1; i>0; i-- ) spaceN = D[i]*spaceN + 1;
    spaceM = (D[i]+1)*spaceN+1;
    spaceP = (2*D[i]-1)*spaceN+1;
    spaceN = D[i]*spaceN+1;
    if( da==0 && ISCONSTANT(A,N) ) {
        z = A[N];
        C[-1] = db;
        for( i=0; i<=db; i++ ) {
            IPCONMUL(z,B,C,N,D,p);
            B += spaceN;
            C += spaceN;
        }
        return;
    }
    if( db==0 && ISCONSTANT(B,N) ) {
        z = B[N];
        C[-1] = da;
        for( i=0; i<=da; i++ ) {
            IPCONMUL(z,A,C,N,D,p);
            A += spaceN;
            C += spaceN;
        }
        return;
    }
    T1 = T; T2 = T+spaceP; T = T2+spaceP;
    dc = da+db;
//if( N==1 && D[0]==2 && M[2]==0 ) z=NEG(M[1],p); else z=0;
//if( z ) { printf("M="); polprint64s(M+1,2,p); printf("\n"); }
    for( k=0; k<=dc; k++ ) {
        T1[0] = -1;
        bot = MAX(0,k-db); top = MIN(k,da);
        for( i=bot; i<=top; i++ )
            if( i==bot ) 
                IPMUL(A+i*spaceN,B+(k-i)*spaceN,T1,N-1,D+1,M+spaceM,T,p);
            else {
                IPMUL(A+i*spaceN,B+(k-i)*spaceN,T2,N-1,D+1,M+spaceM,T,p);
                IPADD(T1,T2,N-1,D+1,p); // T1 += T2
            }
//if( z ) REDUCE(T1,z,p); else
        IPREM(T1,M,N-1,D+1,M+spaceM,T,p);
        COPYAtoB(T1,C+k*spaceN,N,D);
    }
    for( ; dc>=0 && C[dc*spaceN]==-1; dc-- );
    C[-1] = dc; // dc = degree of the product C
    return;
}

void IPREM( LONG *A, LONG *B, int N, int *D, LONG *M, LONG *T, LONG p ) {
// Compute A divided B for A,B in R[x] where B must be monic.
// The remainder r and quotient q are returned in A.
// T = array( TempSpace(N,D) ) is temporary space for products in R
// If N=0, r = 4+5x^2 so dr=deg(r)=2, and q = 7+8x then 
// A = [dr,r,q] = [2,4,0,5,7,8] so deg(q)=1 is not stored in A
// If N==1, D=[2], r=(3+4z)+(5)x^2, q=(6+7z)+(8z)x then
// A = [dr,r,q] = [2,(1,3,4),(-1,*,*),(0,5,*),(1,6,7),(1,0,8)] 
// where * means anything and degree -1 means a 0 coefficient
    int da,db,dr,i,j,k,spaceN,spaceM,spaceP,top,bot,monic; LONG alp, *T1, *T2;
    da = A[0];
    db = B[0];
    A++; B++; // point to coefficients
    if( N==0 ) { A[-1] = poldiv64s(A,B,da,db,p); return; }
//if( N==1 && D[0]==2 && M[2]==0 ) alp=M[1]; else alp=0;
//if( alp ) { printf("M="); polprint64s(M+1,2,p); printf("\n"); }
    if( da<db ) return;
    spaceN = 1;
    for( i=N-1; i>0; i-- ) spaceN = D[i]*spaceN + 1;
    spaceM = (D[0]+1)*spaceN+1;
    spaceP = (2*D[0]-1)*spaceN+1;
    spaceN = D[0]*spaceN + 1;
    T1 = T; T2 = T+spaceP; T = T2+spaceP;
    for( k=da-1; k>=0; k-- ) {
        T1[0] = -1;
        bot = MAX(0,k-(da-db)); top = MIN(db-1,k);
        for( i=bot; i<=top; i++ )
            if( i==bot ) // put in T1
                IPMUL(A+(k-i+db)*spaceN,B+i*spaceN,T1,N-1,D+1,M+spaceM,T,p);
            else {
                IPMUL(A+(k-i+db)*spaceN,B+i*spaceN,T2,N-1,D+1,M+spaceM,T,p);
                IPADD(T1,T2,N-1,D+1,p); // T1 += T2
            }
//if( alp ) REDUCE(T1,alp,p); else
        IPREM(T1,M,N-1,D+1,M+spaceM,T,p);
        IPSUB(A+k*spaceN,T1,N-1,D+1,p);
    }
    for( dr=db-1; dr>=0 && A[dr*spaceN]==-1; dr-- );
    A[-1] = dr; // dr = degree of remainder
    return;
}


int TempSpaceIPINV( int N, int *D ) {
// Space for computing an inverse in R
    int S,I,P;
    if( N==0 ) return 0; // a scalar needs none
    if( N==1 ) return 5*(D[0]+1); // 5 temporaries for polgcdext64s
    S = 4*SpaceMinpoly(N,D); // for r0, r1, t0, t1 in IPINV
    I = TempSpaceIPINV(N-1,D+1); // for recursive inverse in MONIC
    P = TempSpace(N-1,D+1); // for subsequent scalar x in MONIC
    //printf("S=%d I=%d P=%d Total=%d\n",S,I,P,S+MAX(I,P));
    return S+MAX(I,P);
}

int IPINV( LONG *A, int N, int *D, LONG *M, LONG *I, LONG *W, LONG p ) {
// A is in R.
// If A is invertible compute A^(-1) in I and return 0.
// Otherwise return a zero divisor in R and return 1.
// The zero divisor is a factor of one of the minimal polynomials.
    int i,dm,da,spaceN,spaceNm1,spaceM,spaceP,z,d0,d1,dr;
    LONG *r0,*r1,*t0,*t1,*tmp,*T1,*T2;
    if( N==0 ) { printf("IPINV: N must be > 0\n"); exit(1); }
    dm = M[0];
    if( dm!=D[0] ) { printf("M's degree does not match D\n"); exit(1); }
    da = A[0];
    if( da==-1 ) { printf("zero inverse detected\n"); exit(1); }
    if( da>=dm ) { printf("deg(A) must be < deg(M)\n"); exit(1); }
    if( N==1 ) { int dG,dS,dT; LONG *a,*m,*G,*T;
        if( da==0 ) { I[0] = 0; I[1] = modinv64s(A[1],p); return 0; }
        m = W; W += dm+1; polcopy64s(M+1,dm,m);
        a = W; W += dm+1; polcopy64s(A+1,da,a);
        G = W; W += dm+1;
        T = W; W += dm+1;
        // W size must be max(dm+1,da+1) = dm+1
        // So 5 arrays m,a,G,T,W all of size dm+1
        polgcdext64s(m,a,dm,da,G,NULL,T,&dG,&dS,&dT,W,p);
        if( dG==0 ) { polcopy64s(T,dT,I+1); I[0] = dT; return 0; }
        else { polcopy64s(G,dG,I+1); I[0] = dG; return 1; }
    }
    // s_i M + t_i A = r_i
    spaceN = 1;
    for( i=N-1; i>0; i-- ) spaceN = D[i]*spaceN + 1;
    //spaceP = (2*D[0]-1)*spaceN+1;
    spaceNm1 = spaceN;
    spaceM = (D[0]+1)*spaceN+1;
    spaceN = D[0]*spaceN + 1;
    if( da==0 ) { // A is a constant; avoid the EEA
         z = IPINV( A+1, N-1, D+1, M+spaceM, I+1, W, p );
         if( z>0 ) return z; else { I[0] = 0; return 0; }
    }
    // s0 M + t0 A = r0
    r0 = W; W += spaceM; COPYAtoB(M,r0,N,D);  // r0 = M
    t0 = W; W += spaceM; t0[0] = -1;          // t0 = 0
    // s1 M + t1 A = r1
    r1 = W; W += spaceM; COPYAtoB(A,r1,N,D);  // r1 = A
    t1 = W; W += spaceM;
    for( i=0; i<N; i++ ) t1[i] = 0;
    t1[i] = 1;                                // t1 = 1
    M += spaceM;
    d0 = r0[0]; d1 = r1[0];
    while( d1 != -1 ) { LONG *LC;
        z = MONIC( r1, N-1, D+1, M, I+1, W, p );
        if( z>0 ) { I[0] = 0; return z; }
        IPSCAMUL(I+1,t1,N-1,D+1,M,W,p); // t1 = I t1
        IPREM(r0,r1,N-1,D+1,M,W,p);
        dr = r0[0];
        if( dr!=-1 ) { LONG t,*q;
            q = r0+d1*spaceNm1; // point to the quotient in r0
            t = q[0]; // save this value
            q[0] = d0-d1; // set the degree of q
            IPMUL(q,t1,I,N-1,D+1,M,W,p); // I = q t1
            q[0] = t; // restore q[0]
            IPSUB(t0,I,N-1,D+1,p);       // t0 = t0 - q t1
        }
        tmp = r0; r0 = r1; r1 = tmp;
        tmp = t0; t0 = t1; t1 = tmp;
        d0 = d1; d1 = dr;
    }
    if( d0>0 ) { COPYAtoB(r0,I,N,D); return 1; } // zero divisor
    else { COPYAtoB(t0,I,N,D); return 0; } // the inverse
}


int TempSpaceIPGCD( int N, int *D ) {
    int I,P;
    if( N==0 ) return 0; // polgcd64s is inplace
    I = TempSpaceIPINV(N,D);
    P = TempSpace(N,D);
    //printf("Gspace: N=%d I=%d P=%d max=%d\n",N,I,P,MAX(I,P));
    return MAX(I,P);
}

int IPGCD( LONG *A, LONG *B, int N, int *D, LONG *M, LONG *Z, LONG *W, LONG p ) {
    // Input A,B in R[x]
    // If 0 is returned the gcd(A,B) in R[x] is returned in A
    // If >0 is returned a zero divisor in R is returned in Z
    // Both A and B are destroyed note
    // W is an array of working storage = array( TempSpaceIGCD(n,D) ).
    // It needs to be large enough to compute an inverse in R then multiply in R.
    // So max( Ispace(N,D),Mspace(N,D) ) = Ispace(N,D)
    int da,db,dr,z; LONG *r0,*r1,*tmp;
    da = A[0];
    db = B[0];
    if( N==0 ) { A[0] = polgcd64s(A+1,B+1,da,db,p); return 0; }
    r0 = A;
    r1 = B;
    while( db!=-1 ) {
        z = MONIC(r1,N,D,M,Z,W,p);
        if( z>0 ) return z;
        IPREM(r0,r1,N,D,M,W,p); // r0 = r0 div r1
        dr = r0[0]; da = db; db = dr;
        tmp = r0; r0 = r1; r1 = tmp;
    }
    // The gcd is in r0.  It has degree da.  Make it monic.
    z = MONIC(r0,N,D,M,Z,W,p);
    if( z>0 ) return z;
    if( r0!=A ) POLCOPYAtoB(r0,A,N,D);
    return 0;
}


int IPGCDCOF( LONG *A, LONG *B, LONG *CA, LONG *CB, int N, int *D, LONG *M, LONG *Z, LONG *W, LONG p ) {
    int da,db,dg,zd,s,spacerem,i;
    da = A[0];
    db = B[0];
    POLCOPYAtoB( A, CA, N, D );
    POLCOPYAtoB( B, CB, N, D );
    zd = IPGCD( A, B, N, D, M, Z, W, p ); // A = GCD(A,B)
    if( zd ) return zd;
    dg = A[0];
    printf("IPGCDCOF: deg(gcd)=%d\n",dg);
    IPREM(CA,A,N,D,M,W,p);
    if( CA[0]!=-1 ) { printf("GCD did not divide A ?\n"); return -1; }
    IPREM(CB,A,N,D,M,W,p);
    if( CB[0]!=-1 ) { printf("GCD did not divide B ?\n"); return -1; }
    s = SpaceElem(N,D);
    CA[0] = da-dg; // deg(A/G)
    CB[0] = db-dg; // deg(B/G)
    spacerem = dg*s;
    for( i=1; i<=(da-dg+1)*s; i++ ) CA[i] = CA[i+spacerem];
    for( i=1; i<=(db-dg+1)*s; i++ ) CB[i] = CB[i+spacerem];
    return 0;
}
    

void IPremovedegreesRING( LONG *A, LONG *B, int N, int *D ) {
    int d,i,j,s,t;
    if( N<1 ) { printf("IPremovedegreesRING: N must be > 0"); exit(1); }
    d = A[0];
    A++;
    if( N==1 ) { 
        for( i=0; i<=d; i++ ) B[i] = A[i]; 
        for( i=d+1; i<D[0]; i++ ) B[i] = 0;
        return;
    }
    s = SpaceElem(N-1,D+1);
    for( t=1,i=1; i<N; i++ ) t = t*D[i]; // dense space
    for( i=0; i<=d; i++, A+=s,B+=t ) IPremovedegreesRING( A, B, N-1, D+1 );
    for( i=d+1; i<D[0]; i++,B+=t ) for( j=0; j<t; j++ ) B[j] = 0;
    return;
}

void IPremovedegreesPOLY( LONG *A, LONG *B, int N, int *D ) {
// Makes the coefficients of A in R[x] dense 
// E.g. if A=[1,1,2,3,0,2,4,5,6] = (2+3*y)+(4+5*y+6*y^2)*x and N=1, D=[3]
// then B = [2,3,0,4,5,6].  Note, B may equal A.
    int d,i,s,t;
    d = A[0]; A++;
    if( N==0 ) { 
        for( i=0; i<=d; i++ ) B[i] = A[i];
        return;
    }
    s = SpaceElem(N,D); // space for R
    for( t=1,i=0; i<N; i++ ) t = t*D[i]; // dense space
    for( i=0; i<=d; i++ ) IPremovedegreesRING( A+i*s, B+i*t, N, D );
    return;
}

void IPinsertdegreesRING( LONG *A, LONG *B, int N, int *D ) {
    int d,i,s,t;
    if( N<1 ) { printf("IPinsertdegreesRING: N must be > 0"); exit(1); }
    if( N==1 ) { 
        for( d=D[0]-1; d>=0 && A[d]==0; d-- );
        B++;
        for( i=d; i>=0; i-- ) B[i] = A[i]; // go backwards!
        B[-1] = d;
        return;
    }
    s = SpaceElem(N-1,D+1);
    for( t=1,i=1; i<N; i++ ) t = t*D[i]; // dense space
    d = D[0]-1; B++;
    for( i=d; i>=0; i-- ) IPinsertdegreesRING( A+i*t, B+i*s, N-1, D+1 ); // go backwards!
    while( d>=0 && B[d*s]==-1 ) d--;
    B[-1] = d;
    return;
}

void IPinsertdegreesPOLY( LONG *A, int d, LONG *B, int N, int *D ) {
// Makes the coefficients of A in R[x] recursive dense
// E.g. if A=[2,3,0,4,5,6], d=1, N=1, D=[3] so that A = (2+3*y)+(4+5*y+6*y^2)*x
// then B = [1,1,2,3,0,2,4,5,6].  Note, B may equal A.
    int i,s,t;
    if( N==0 ) {
        while( d>=0 && A[d]==0 ) d--;
        B[0] = d; B++;
        for( i=0; i<=d; i++ ) B[i] = A[i];
        return;
    }
    s = SpaceElem(N,D); // space for R
    for( t=1,i=0; i<N; i++ ) t = t*D[i]; // dense space
    B++;
    for( i=d; i>=0; i-- ) IPinsertdegreesRING( A+i*t, B+i*s, N, D ); // go backwards!
    while( d>=0 && B[d*s]==-1 ) d--;
    B[-1] = d;
    return;
}


void matvecmul64s( LONG *A, int n, LONG *x, LONG *y, LONG p );
LONG matinv64s( LONG * A, int n, LONG * B, LONG p );
void matprint64s( LONG *A, int n );

LONG IPmakephi( 
   LONG *gamma, int N, int *D, LONG *M,
   LONG *phi, LONG *I, LONG *m, // outputs
   LONG p ) {
// There is a linear relation between { gamma^i : 0 <= i <= d }
// Build the multiplication matrix phi where col(i,phi) = gamma^(i-1)
// Compute phi^(-1) in I and use I to compute the minimal polynomial in m
// Let d = [R:Zp] = deg(R).  Then phi is dxd and m has degree d.
// Output det(phi).  If det(phi)=0 then the mapping does not exist.
    int d,i,k,spaceN,spaceT; 
    LONG det,*g1,*g2,*g3,*T,*A;
//clock_t st,tt,pt,mt;
//tt = clock();
    for( d=1,i=0; i<N; i++ ) d = d*D[i]; // d = degree(R)
//printf("IPmakePhi: N=%d  d=%d\n",N,d);
//printf("gamma = "); printpol( gamma, N-1, D+1, "xyz" ); printf("\n");
    for( phi[0]=1,i=1; i<d; i++ ) phi[i*d] = 0; // col(1,phi) = 1
    IPremovedegreesRING(gamma,m,N,D);
    for( i=0; i<d; i++ ) phi[i*d+1] = m[i]; // col(2,phi) = gamma
    spaceN = SpaceElem(N,D);
    g2 = array(spaceN+1); if( debug ) g2[spaceN] = 2222;
    g3 = array(spaceN+1); if( debug ) g3[spaceN] = 3333;
    COPYAtoB(gamma,g2,N,D);
    if( debug && g2[spaceN]!=2222 ) printf("g2 mem bad\n");
    spaceT = TempSpace(N,D);
    T = array(spaceT+1); if( debug ) T[spaceT] = 1111;
//mt = clock();
    for( k=2; k<d; k++ ) {
        IPELTMUL(gamma,g2,g3,N,D,M,T,p); // g3 = gamma g2 = gamma^k
        if( debug && T[spaceT] !=1111 ) printf(" IPmakephi T overflow\n");
        IPremovedegreesRING(g3,m,N,D); // compress to dense
        for( i=0; i<d; i++ ) phi[i*d+k] = m[i]; // col(k,phi) = gamma^k
        g1 = g2; g2 = g3; g3 = g1;
    }
//pt = clock();
//mt = pt-mt;
    A = array(d*d+1); A[d*d] = 4444; // temporary storage
    veccopy64s(phi,d*d,A);
    det = matinv64s(A,d,I,p); // I = A^(-1)
//st = clock();
//pt = st-pt;
    if( debug && A[d*d]!=4444 ) printf("A mem bad\n");
    if( det==0 ) {
        free(g2); free(g3); free(T); free(A);
        return 0; // matrix is not invertible so stop
    }
    IPELTMUL(gamma,g2,g3,N,D,M,T,p); // g3 = gamma^d
    IPremovedegreesRING(g3,g3,N,D); // m = gamma^d
    m[0] = d;
    m++;
    matvecmul64s( I, d, g3, m, p ); // m = I g3 = phi^(-1) gamma^d

    for( i=0; i<d; i++ ) m[i] = NEG(m[i],p); // m = -m
    m[d] = 1;
    if( debug && T[spaceT]!=1111 ) printf("T mem bad\n");
    if( debug && g2[spaceN]!=2222 ) printf("g2 mem bad\n");
    if( debug && g3[spaceN]!=3333 ) printf("g3 mem bad\n");
    free(g2); free(g3); free(T); free(A);
//tt = clock()-tt;
//printf("phi:  tot=%.3fms  phi=%.3fms  mul=%.3fms\n",
//    tt/1000., pt/1000., mt/1000. );
    return det;
}


int IPphiGCD( LONG *A, LONG *B, int N, int *D, LONG *M, 
            LONG *phi, LONG *Iphi, LONG *minpoly,
            LONG *Z, LONG *W, LONG p ) {
// Same as PHIGCD but takes phi and it's inverse and m as inputs
    int i,d,da,db,dg,DD[1];

    if( N==1 ) return IPGCD(A,B,N,D,M,Z,W,p);
    for( d=1,i=0; i<N; i++ ) d = d*D[i];

    // apply matrix Iphi to coefficients of A and B in R
    da = A[0]; IPremovedegreesPOLY(A,A,N,D);
    db = B[0]; IPremovedegreesPOLY(B,B,N,D);
    LONG *u = W; // use W for the vector u
    for( i=0; i<=da; i++ ) {
        matvecmul64s(Iphi,d,A+i*d,u,p); // u = phi coeff(A,x,i)
        veccopy64s(u,d,A+i*d); }
    for( i=0; i<=db; i++ ) {
        matvecmul64s(Iphi,d,B+i*d,u,p); // u = phi coeff(B,x,i)
        veccopy64s(u,d,B+i*d); }

    DD[0] = d;
    IPinsertdegreesPOLY(A,da,A,1,DD);
    IPinsertdegreesPOLY(B,db,B,1,DD);

    i = IPGCD(A,B,1,DD,minpoly,Z,W,p);

    if( i>0 ) { // zero divisor
        IPremovedegreesRING(Z,Z,1,DD);
        matvecmul64s(phi,d,Z,u,p);
        IPinsertdegreesRING(u,Z,N,D);
        return 1;
    }
    dg = A[0];

    IPremovedegreesPOLY(A,A,1,DD);
    for( i=0; i<=dg; i++ ) {
        matvecmul64s(phi,d,A+i*d,u,p);
        veccopy64s(u,d,A+i*d); }
    IPinsertdegreesPOLY(A,dg,A,N,D);
    return 0;
}


int IPphiGCDCOF( LONG *A, LONG *B, LONG *CA, LONG *CB, int N, int *D, LONG *M, 
            LONG *phi, LONG *Iphi, LONG *minpoly,
            LONG *Z, LONG *W, LONG p ) {
// Same as PHIGCD but takes phi and it's inverse and m as inputs
    int i,d,da,db,dg,s,spacerem,DD[1];
    clock_t st,tt,gcdt,phit,remt;

    if( N==1 ) return IPGCD(A,B,N,D,M,Z,W,p);
    for( d=1,i=0; i<N; i++ ) d = d*D[i];

tt = st = clock();
    // apply matrix Iphi to coefficients of A and B in R
    da = A[0]; IPremovedegreesPOLY(A,A,N,D);
    db = B[0]; IPremovedegreesPOLY(B,B,N,D);
    LONG *u = W; // use W for the vector u
    for( i=0; i<=da; i++ ) {
        matvecmul64s(Iphi,d,A+i*d,u,p); // u = phi coeff(A,x,i)
        veccopy64s(u,d,A+i*d); }
    for( i=0; i<=db; i++ ) {
        matvecmul64s(Iphi,d,B+i*d,u,p); // u = phi coeff(B,x,i)
        veccopy64s(u,d,B+i*d); }
phit = clock()-st;

    DD[0] = d;
    IPinsertdegreesPOLY(A,da,A,1,DD);
    IPinsertdegreesPOLY(B,db,B,1,DD);

    da = A[0];
    db = B[0];
    POLCOPYAtoB( A, CA, 1, DD );
    POLCOPYAtoB( B, CB, 1, DD );

st = clock();
    i = IPGCD(A,B,1,DD,minpoly,Z,W,p);
gcdt = clock()-st;

    if( i>0 ) { // zero divisor
        IPremovedegreesRING(Z,Z,1,DD);
        matvecmul64s(phi,d,Z,u,p);
        IPinsertdegreesRING(u,Z,N,D);
        return 1;
    }
    dg = A[0];

st = clock();
    IPREM(CA,A,1,DD,minpoly,W,p);
    if( CA[0]!=-1 ) { printf("GCD did not divide A ?\n"); return -1; }
    IPREM(CB,A,1,DD,minpoly,W,p);
    if( CB[0]!=-1 ) { printf("GCD did not divide B ?\n"); return -1; }
remt = clock()-st;

    s = SpaceElem(1,DD);
    CA[0] = da-dg; // deg(A/G)
    CB[0] = db-dg; // deg(B/G)
    spacerem = dg*s;
    for( i=1; i<=(da-dg+1)*s; i++ ) CA[i] = CA[i+spacerem];
    for( i=1; i<=(db-dg+1)*s; i++ ) CB[i] = CB[i+spacerem];

    IPremovedegreesPOLY(A,A,1,DD);
    IPremovedegreesPOLY(CA,CA,1,DD);
    IPremovedegreesPOLY(CB,CB,1,DD);

st = clock();
    for( i=0; i<=dg; i++ ) {
        matvecmul64s(phi,d,A+i*d,u,p);
        veccopy64s(u,d,A+i*d); }
    IPinsertdegreesPOLY(A,dg,A,N,D);
    for( i=0; i<=da-dg; i++ ) {
        matvecmul64s(phi,d,CA+i*d,u,p);
        veccopy64s(u,d,CA+i*d); }
    IPinsertdegreesPOLY(CA,da-dg,CA,N,D);
    for( i=0; i<=db-dg; i++ ) {
        matvecmul64s(phi,d,CB+i*d,u,p);
        veccopy64s(u,d,CB+i*d); }
    IPinsertdegreesPOLY(CB,db-dg,CB,N,D);
phit += clock()-st;
tt = clock()-tt;

printf("PHIGCDCOF:  time=%.3fs  phi=%.3fs(%.1f%%)  gcd=%.3fs(%.1f%%)  rem=%.3fs(%.1f%%)\n",
    tt/1000000., phit/1000000., 100.0*phit/tt, gcdt/1000000., 100.0*gcdt/tt, remt/1000000., 100.0*remt/tt );
    return 0;
}


int IPphiDIV2( LONG *A, LONG *B, LONG *G, int N, int *D, LONG *M, 
            LONG *phi, LONG *Iphi, LONG *minpoly,
            LONG *Z, LONG *W, LONG p ) {

// Test if G|A and G|B in R[x] using phi

    int i,d,da,db,dg,s,spacerem,DD[1];
    clock_t st,tt,gcdt,phit,remt;

    for( d=1,i=0; i<N; i++ ) d = d*D[i];

tt = st = clock();
    // apply matrix Iphi to coefficients of A and B in R
    da = A[0]; IPremovedegreesPOLY(A,A,N,D);
    db = B[0]; IPremovedegreesPOLY(B,B,N,D);
    dg = G[0]; IPremovedegreesPOLY(G,G,N,D);

    LONG *u = W; // use W for the vector u
    for( i=0; i<=da; i++ ) {
        matvecmul64s(Iphi,d,A+i*d,u,p); // u = phi coeff(A,x,i)
        veccopy64s(u,d,A+i*d); }
    for( i=0; i<=db; i++ ) {
        matvecmul64s(Iphi,d,B+i*d,u,p); // u = phi coeff(B,x,i)
        veccopy64s(u,d,B+i*d); }
    for( i=0; i<=dg; i++ ) {
        matvecmul64s(Iphi,d,G+i*d,u,p); // u = phi coeff(G,x,i)
        veccopy64s(u,d,G+i*d); }
phit = clock()-st;

    DD[0] = d;
    IPinsertdegreesPOLY(A,da,A,1,DD);
    IPinsertdegreesPOLY(B,db,B,1,DD);
    IPinsertdegreesPOLY(G,dg,G,1,DD);

    da = A[0];
    db = B[0];
    dg = G[0];

st = clock();
    IPREM(A,G,1,DD,minpoly,W,p);
    if( A[0]!=-1 ) { printf("GCD did not divide A ?\n"); return 0; } 
    //else printf("G|A\n");
    IPREM(B,G,1,DD,minpoly,W,p);
    if( B[0]!=-1 ) { printf("GCD did not divide B ?\n"); return 0; }
    //else printf("G|B\n");
remt = clock()-st;

    s = SpaceElem(1,DD);
    A[0] = da-dg; // deg(A/G)
    B[0] = db-dg; // deg(B/G)
    spacerem = dg*s;
    // IPREM outputs A rem G (of degree dg-1) followed by and A quo G in A
    // Move the quotient (of degree da-dg) to the front of A.  Same for B.
    for( i=1; i<=(da-dg+1)*s; i++ ) A[i] = A[i+spacerem];
    for( i=1; i<=(db-dg+1)*s; i++ ) B[i] = B[i+spacerem];

    IPremovedegreesPOLY(A,A,1,DD);
    IPremovedegreesPOLY(B,B,1,DD);

st = clock();
    for( i=0; i<=da-dg; i++ ) {
        matvecmul64s(phi,d,A+i*d,u,p);
        veccopy64s(u,d,A+i*d); }
    for( i=0; i<=db-dg; i++ ) {
        matvecmul64s(phi,d,B+i*d,u,p);
        veccopy64s(u,d,B+i*d); }

    IPinsertdegreesPOLY(A,da-dg,A,N,D);
    IPinsertdegreesPOLY(B,db-dg,B,N,D);
phit += clock()-st;
tt = clock()-tt;

printf("PHICOF2:  time=%.3fs  phi=%.3fs(%.1f%%)  rem=%.3fs(%.1f%%)\n",
    tt/1000000., phit/1000000., 100.0*phit/tt, remt/1000000., 100.0*remt/tt );
    return 1;
}


int PHIGCD( LONG *A, LONG *B, int N, int *D, LONG *M, 
            LONG *gamma,
            LONG *Z, LONG *W, LONG p ) {
// Same as IPGCD but use primitive element gamma in R to map R into Zp[z]/m(z)
    int i,d,da,db,dg,DD[1];
    LONG det,*phi,*Iphi,*minpoly;
    clock_t st,tott,gcdt,phit,rt;

tott = clock();

    if( N==1 ) return IPGCD(A,B,N,D,M,Z,W,p);
    for( d=1,i=0; i<N; i++ ) d = d*D[i];

printf("PHIGCD: da=%lld  db=%lld  d=%d\n",A[0],B[0],d);
//printf("A="); printpol(A, N, D, "xyzuvw" ); printf("\n");
//printf("B="); printpol(B, N, D, "xyzuvw" ); printf("\n");

    phi = array(d*d+1); phi[d*d] = 1111;
    Iphi = array(d*d+1); Iphi[d*d] = 2222; 
    minpoly = array(d+3); minpoly[d+2] = 3333;

    det = IPmakephi( gamma, N, D, M, phi, Iphi, minpoly, p );
    if( det==0 ) {printf("gamma is no good\n"); return -1; }
//printf("phi:\n"); matprint64s( phi, d ); printf("\n");
//printf("invphi:\n"); matprint64s( Iphi, d ); printf("\n");
//printf("M="); printpol(minpoly,0,D,"zxy"); printf("\n");


if( debug && phi[d*d] != 1111 ) printf("phi mem bad\n");
if( debug && Iphi[d*d] != 2222 ) printf("Iphi mem bad\n");
if( debug && minpoly[d+2] != 3333 ) printf("minpoly mem bad\n");

st = clock(); phit = st-tott;

    i = IPphiGCD( A, B, N, D, M, phi, Iphi, minpoly, Z, W, p );

rt = clock(); gcdt = rt-st; tott = rt-tott;
printf("PHIGCD:  time=%.3fs  gcd=%.3fs  phi=%.3fs(%.1f%%)\n",
    tott/1000000.0, gcdt/1000000.0, phit/1000000.0, 100.0*phit/tott );

    free(phi); free(Iphi); free(minpoly);
    return i;
}


int PHIGCDCOF( LONG *A, LONG *B, LONG *CA, LONG *CB,
               int N, int *D, LONG *M, LONG *gamma,
               LONG *Z, LONG *W, LONG p ) {
// Computes G = GCD(A,B) in A and CA = A/G and CB = B/G
// Same as IPGCDCOF but use primitive element gamma in R to map R into Zp[z]/m(z)
    int i,d,dg;
    LONG det,*phi,*Iphi,*minpoly;
    clock_t st,tott,gcdt,phit,rt;

tott = clock();

    if( N==1 ) return IPGCDCOF(A,B,CA,CB,N,D,M,Z,W,p);
    for( d=1,i=0; i<N; i++ ) d = d*D[i];

printf("PHIGCDCOF: da=%lld  db=%lld  d=%d\n",A[0],B[0],d);
//printf("A="); printpol(A, N, D, "xyzuvw" ); printf("\n");
//printf("B="); printpol(B, N, D, "xyzuvw" ); printf("\n");

    phi = array(d*d+1); phi[d*d] = 1111;
    Iphi = array(d*d+1); Iphi[d*d] = 2222; 
    minpoly = array(d+3); minpoly[d+2] = 3333;

    det = IPmakephi( gamma, N, D, M, phi, Iphi, minpoly, p );
    if( det==0 ) return -1; // gamma is no good

if( debug && phi[d*d] != 1111 ) printf("phi mem bad\n");
if( debug && Iphi[d*d] != 2222 ) printf("Iphi mem bad\n");
if( debug && minpoly[d+2] != 3333 ) printf("minpoly mem bad\n");

st = clock(); phit = st-tott;

    i = IPphiGCDCOF( A, B, CA, CB, N, D, M, phi, Iphi, minpoly, Z, W, p );

rt = clock(); gcdt = rt-st; tott = rt-tott;
printf("PHIGCDCOF:  time=%.3fs  gcd=%.3fs  phi=%.3fs(%.1f%%)\n",
    tott/1000000.0, gcdt/1000000.0, phit/1000000.0, 100.0*phit/tott );

    free(phi); free(Iphi); free(minpoly);
    return i;
}


int PHICOF2( LONG *A, LONG *B, LONG *G,
               int N, int *D, LONG *M, LONG *gamma,
               LONG *Z, LONG *W, LONG p ) {
// Computes A = A/G and B = B/G.  Returns 1 if both divide, 0 otherwise.
// Same as IPGCDCOF but use primitive element gamma in R to map R into Zp[z]/m(z)
    int i,d;
    LONG det,*phi,*Iphi,*minpoly;
    clock_t st,tott,divt,phit,rt;

tott = clock();

    for( d=1,i=0; i<N; i++ ) d = d*D[i];

printf("PHICOF2: da=%lld  db=%lld  dg=%lld  d=%d\n",A[0],B[0],G[0],d);

    phi = array(d*d+1); phi[d*d] = 1111;
    Iphi = array(d*d+1); Iphi[d*d] = 2222; 
    minpoly = array(d+3); minpoly[d+2] = 3333;

    det = IPmakephi( gamma, N, D, M, phi, Iphi, minpoly, p );
    if( det==0 ) return -1; // gamma is no good

if( debug && phi[d*d] != 1111 ) printf("phi mem bad\n");
if( debug && Iphi[d*d] != 2222 ) printf("Iphi mem bad\n");
if( debug && minpoly[d+2] != 3333 ) printf("minpoly mem bad\n");

st = clock(); phit = st-tott;

    i = IPphiDIV2( A, B, G, N, D, M, phi, Iphi, minpoly, Z, W, p );

rt = clock(); divt = rt-st; tott = rt-tott;
printf("PHICOF2:  time=%.3fs  div=%.3fs  phi=%.3fs(%.1f%%)\n",
    tott/1000000.0, divt/1000000.0, phit/1000000.0, 100.0*phit/tott );

    free(phi); free(Iphi); free(minpoly);
    return i;
}






/*********************************************************************/
/*  These routines work with the Maple DAG directly                  */
/*********************************************************************/
// Copyright, Michael Monagan, April 2022, May 2022

#define NAME 8
#define TABLEREF 10
#define PROD 14
#define SUM 16
#define POLY 17
#define SEQ 31
#define LIST 32
#define SET 33
#define INTNEG 1
#define INTPOS 2


int ID( LONG *a ) {
    int id;
    id = (a[0] >> 32) & 63;
    return id;
}

LONG LENGTH( LONG *a ) {
    LONG len;
    len = a[0] & ((1LL << 32)-1);
    return len;
}


void recunpackn( LONG *L, int N, int *D, LONG *A, LONG p )
{   // L encodes a polynonial in F[x], F a number field mod p
    // D is the degree of N minimal polynomials e.g. for <y^3+z+5,z^2+7>, D = [3,2]
    // Here each element in F needs D[0]*(D[1]+1)+1 = 3*(2+1)+1 = 10 words
    // Upack L into the array A using a dense recursive encoding with degrees in A
    // E.g. for f = (5+8z+7*y^2)x^1 we have L = [0,[[5,8],0,[7]]] hence
    //   A = [1, (-1,*,*,*,*,*,*,*,*,*), (2,(1,5,8),(-1,*,*),(0,7,*))]
    // where * means don't care and degree=-1 means a 0 coefficient.
    int i,d,sn;
    LONG c = (LONG) L;
    if( c&1 ) {
        if( N==-1 ) { c = c>>1; A[0] = mod64s(c,p); return; }
        else if( c==1 ) { A[0] = -1; return; }
        else { printf("integer unexpected\n"); exit(1); }
    }
    sn = 1; // sn = space for an element of F
    for( i=N-1; i>=0; i-- ) sn = D[i]*sn+1;
    if( ID(L)==LIST ) L = (LONG *) L[1]; // to SEQ
    else { printf("recunpack: list expected\n"); exit(1); }
    d = LENGTH(L)-2;
    A++;
    L++;
    for( i=0; i<=d; i++ ) recunpackn((LONG *) L[i],N-1,D+1,A+i*sn,p);
    while( d>=0 && A[d*sn]==-1 ) d--; // check the degree!
    A[-1] = d;
    return;
}

void UnpackMinPolys( LONG *L, int N, int *D, LONG *A, LONG p ) {
    int i,k,sn; LONG c;
    if( ID(L)!=LIST ) { printf("minpolys: list expected\n"); exit(1); }
    L = (LONG *) L[1]; // point to SEQ
    if( LENGTH(L)!=N+1 ) { printf("minpolys: bad length\n"); exit(1); }
    for( k=1; k<LENGTH(L); k++ ) {
        LONG *M = (LONG *) L[k];
        if( ID(M)!=LIST ) { printf("minpolys: list expected\n"); exit(1); }
        recunpackn( M, N-1, D+1, A, p );
        for( sn=1,i=N-1; i>=0; i-- ) // compute sn = space for M
            sn = i==0 ? (D[i]+1)*sn+1 : D[i]*sn+1; // +1 for degree
        A += sn; D++; N--; // to next minpoly
    }
    return;
}