/* Bivariate Hensel lift - cubic cost */
/* Garrett Paluck, 2021 */
//V33

#include<stdio.h>
#include<stdlib.h>
#include<time.h>
#include<math.h>
#include <sys/time.h>

#define DEBUG 0

#define M_INT long long int
#define LONG long long int
#define ULONG unsigned long long int
#define ULONG unsigned long long int
#define UINT32 unsigned int
#define UINT64 unsigned long long
typedef enum{false, true} bool;

#define WORDSIZE 64
typedef unsigned long M_UINT;

/* [x0:x1] := a*b */
#define MUL211(x0,x1,a,b) __asm__(\
"       mulq    %3      \n\t"   \
: "=a"(x0), "=d"(x1) : "0"(a), "r"(b) : "cc")

/* x0 := [x0:x1] / a */
/* x1 := [x0:x1] % a */
#define DIV21H(x0,x1,a) __asm__(\
"       divq    %2      \n\t"   \
: "=a"(x0), "=d"(x1) : "r"(a), "0"(x0), "1"(x1) : "cc")

typedef struct {
	UINT64 s;	/* shift */
	UINT64 v;	/* reciprocal of d */
	UINT64 d0;	/* divisor shifted up */
	UINT64 d1;
} recint;

recint recip1(UINT64  p);
ULONG mulrec64(ULONG  a, ULONG  b, recint  v);

        /* z += a1:a0 */
        #define zadd(z,a0,a1) __asm__(\
        "       addq    %4, %0  \n\t" \
        "       adcq    %5, %1  \n\t" \
                : "=&r"(z[0]), "=&r"(z[1]) : "0"(z[0]), "1"(z[1]), "r"(a0), "r"(a1))

/* z = a*b */
        #define zmul(z,a,b) __asm__(\
        "       mulq    %%rdx   \n\t" \
                : "=a"(z[0]), "=d"(z[1]) : "a"(a), "d"(b))

/* z += a*b */
        #define zfma(z,a,b) do {        \
        unsigned long u,v;              \
        __asm__(                        \
        "       mulq    %%rdx           \n\t" \
        "       addq    %%rax, %0       \n\t" \
        "       adcq    %%rdx, %1       \n\t" \
                : "=&r"(z[0]), "=&r"(z[1]), "=a"(u), "=d"(v) : "0"(z[0]), "1"(z[1]), "a"(a), "d"(b));\
        } while (0)

/* z -= a*b */
        #define zfms(z,a,b) do {        \
        unsigned long u,v;              \
        __asm__(                        \
        "       mulq    %%rdx           \n\t" \
        "       subq    %%rax, %0       \n\t" \
        "       sbbq    %%rdx, %1       \n\t" \
                : "=&r"(z[0]), "=&r"(z[1]), "=a"(u), "=d"(v) : "0"(z[0]), "1"(z[1]), "a"(a), "d"(b));\
        } while (0)

/* z = z % p safe */
        #define zmod(z,p) __asm__(\
        "       divq    %4      \n\t" \
        "       xorq    %0, %0  \n\t" \
                : "=a"(z[1]), "=d"(z[0]) : "a"(z[0]), "d"(z[1] < p ? z[1] : z[1] % p), "r"(p))

ULONG seed, mult;
LONG rand64s(LONG p);

int max(int a, int b);
LONG * array(int n) { return (LONG *) malloc(n*sizeof(LONG)); }
int * arrayI(int n) { return (int *) malloc(n*sizeof(int)); }

inline int max32s(int a, int b) { if(a>b) return a; else return b; }
inline int min32s(int a, int b) { if(a<b) return a; else return b; }
//inline LONG rand64s(LONG p) { LONG x,y; seed = mult*seed; x = seed >> 32; seed = mult*seed; y = seed >> 32; x = (x<<31) | y; x = x % p; return(x);}
inline LONG add64s(LONG a, LONG b, LONG p) { LONG t; t = (a-p)+b; t += (t>>63) & p; return t; }
inline void add264s(LONG a, LONG b, LONG *c, LONG p) {LONG t; t = (a-p)+b; t += (t>>63) & p; *c=t;}
inline LONG sub64s(LONG a, LONG b, LONG p) { LONG t; t = a-b; t += (t>>63) & p; return t; }
inline LONG mul64s(LONG a, LONG b, LONG p) {
        LONG q, r;
        __asm__ __volatile__(           \
        "       mulq    %%rdx           \n\t" \
        "       divq    %4              \n\t" \
        : "=a"(q), "=d"(r) : "0"(a), "1"(b), "rm"(p));
        return r;
}
LONG modinv64s( LONG c, LONG p );

LONG log2N(LONG n) {

   LONG i,j;

   i=0;j=1;
   while(j<n){
      i++; j *= 2;
   }
   return i;

}

LONG ceil2(LONG n){
  if (n%2 == 0){ return n/2;}
  else if (n==1){return 0;}
  else {return (n/2)+1;}
}

void polprint64s( LONG *A, int d );
int poldeg64s(LONG *A, int d);
void polcopy64s( LONG *A, int d, LONG *B );
int poladd64s(LONG *a, LONG *b, LONG *c, int da, int db, LONG p);
void poladdFast(LONG *A,LONG *B, LONG *C,int dx, LONG p);
int polsub64s(LONG *a, LONG *b, LONG *c, int da, int db, LONG p);
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 );
LONG poleval64s(LONG *a, int d, LONG x, LONG p);
int polbivareval64s( LONG *A, LONG *X, int na, LONG *W, LONG *B, int dx, int dy, LONG p );
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 );
int poladdsca64s(LONG *a, LONG *b, LONG *c, int da, int db, LONG alpha, LONG p, recint P);
int polsubsca64s(LONG *a, LONG *b, LONG *c, int da, int db, LONG alpha, LONG p, recint P);


//Custom Functions

//fills an array with 0s
void cleanArray64s(LONG *A, int len);
//prints bivariate Polynomials in Maple orientation
 void polprintbivar64s( LONG *A, int dx, int dy, char * x, char * y);
//Expands a polynomial with dy factors (y-alpha)
void ExpandOutBivar(LONG *A,int dx, int dy, LONG alpha, LONG p);
//Compare solution output by the Cubic function to the actual answer
void ComparePolyFactorsBivar(LONG *A,LONG *F,int dx,int dz,int n,LONG alpha,LONG p);

LONG DotProduct(LONG *A, LONG *B, int d, LONG p){

   int i,j,k;
   ULONG z[2];

   z[0] = 0; z[1] = 0;

   i=0;

   while(i<=d-3){
      zfma(z,A[i],B[i]);i++;
      zfma(z,A[i],B[i]);i++;
      zfma(z,A[i],B[i]);i++;
   }

   while(i<=d){
      zfma(z,A[i],B[i]);i++;
   }
   zmod(z,p);
   return z[0];

}

LONG DotProductP(LONG *A, LONG *B, int d, LONG p){

   int i,j,k;
   ULONG z[2];

   z[0] = 0; z[1] = 0;

   i=0;

   while(i<=d-2){
      zfma(z,A[i],B[i]);i++;//if(z[1]>=p){z[1] -=p;}
      zfma(z,A[i],B[i]);i++;if(z[1]>=p){z[1] -=p;}
   }

   while(i<=d){
      zfma(z,A[i],B[i]);i++;if(z[1]>=p){z[1] -=p;}
   }
   zmod(z,p);
   return z[0];

}

ULONG fastPolyMult(LONG *poly1,LONG *poly2, LONG *T, int deg1, int deg2, LONG p){

    //local variables
    int k,m,min,max;
    ULONG z[2];

    //clean T
    //cleanArray64s(T,deg1+deg2);

    for(k=0;k<=deg1+deg2;k++){
       min = max32s(0,k-deg2);
       max = min32s(k,deg1);
       z[0] = 0; z[1] = 0;
       m = min;
       while(m<=max-3){
          zfma(z,poly1[m],poly2[k-m]);m++;
          zfma(z,poly1[m],poly2[k-m]);m++;
          zfma(z,poly1[m],poly2[k-m]);m++;
       }
       while(m<=max){
          zfma(z,poly1[m],poly2[k-m]);m++;
       }
       zmod(z,p);
       T[k] = add64s(z[0],T[k],p);
     }
}

//perform the Extended Euclidean Algorithm
int MultiEEA(LONG *U, int d, int n, LONG *M, LONG *W, LONG p){
//size of U = n*(d+1)
//size of M = ...
//size of W = 2*((n-1)*d+1) + d+1

   //local variables
   int i,j,du,dm,ds,dt,dg,dmCalc,mPos,point;
   LONG *u,*m,*mCalc,*g,*Wptr,t,alpha2;

   //clean array
   cleanArray64s(W,4*(n*d+1) + (d+1));

   //set initial array pointers
   mCalc = W;
   g = mCalc + (n-1)*d+1;

   //define initial position in the M array
   mPos = 0;
   for(i=n-1;i>1;i--){
      mPos = mPos + i*d+1;
   }

   mCalc[0] = 1;
   dmCalc = 0;

   //Main loop
   u = U + (n-1)*(d+1);
   for (i=n;i>=2;i--){

      //calculate mcalc*U_i
      dmCalc = polmul64s(mCalc,u,mCalc,dmCalc,d,p);
      //add new mcalc to M
      polcopy64s(mCalc,dmCalc,M+mPos);
      mPos = mPos - ((n-i+2)*d + 1);
      //update u pointer
      u -= d+1;
   }


   //make sure gcd(u_i,M_i) = 1 for all M
   u = g + (n-1)*d+1;
   Wptr = u + d+1;
   dmCalc = (n-1)*d;
   mPos = 0;

   for(i=0;i<n;i++){

      //move data to u,mCalc
      polcopy64s(U+i*(d+1),d,u);
      polcopy64s(M+mPos,dmCalc,mCalc);
      polgcdext64s(u,mCalc,d,dmCalc,g,0,0,&dg,&ds,&dt,Wptr,p);

      //fail state
      if(dg>0) return -1;


      dmCalc -= d;
      mPos += (n-i-1)*d+1;
   }

   return 0;

}

void addPol(LONG *A, int dz, LONG p);

void altDiv2(LONG *A,int dx, int dz, LONG alpha, LONG *W, LONG p);

void generateGcdexS(LONG *U, int n, int du, LONG *M, int dx, LONG *S, int *sDeg, LONG *W, LONG p){

    //local variables
    int i,j,k,du2,dm,dg,ds,wTemp,dt,mPos;
    LONG *uTemp,*mTemp,*g,*wPtr;

    cleanArray64s(W,4*n*(du+1));

    uTemp = W;
    mTemp = uTemp + du+1;
    g = mTemp + (n-1)*du+1;
    wPtr = g + n*du+1;

    mPos = 0;

    //Calculate the sigma term and subtract it from C
    for (i=1;i<n;i++){

       //extract M_i,U_i
       polcopy64s(U+(i-1)*(du+1),du,uTemp);
       du2 = poldeg64s(uTemp,du);
       polcopy64s(M+mPos,(n-i)*du,mTemp);
       dm = poldeg64s(mTemp,(n-i)*du);
       mPos = mPos + (n-i)*du+1;

       //perform EEA to solve s*U[i] + t*M[i] = 1 for s in Z[x]
       polgcdext64s(mTemp,uTemp,dm,du2,g,S+(i-1)*(dx+1),0,&dg,&ds,&dt,wPtr,p);
       sDeg[i-1] = ds;

    }

}

void Dionvar(LONG *U, int n, int d, LONG *c, int dx, LONG *M, LONG *S, int *sDeg, LONG *W, LONG *Sigmas, LONG p, LONG *numMuls){
//Enter list of polynomials U of size (n*(d+1)), solution C of size at most (d+1), and polynomial M <- U[2]*U[3]...U[n]

//local variables
    int i,j,k,du,dm,dg,ds,dt,mPos,dck,drem,point;
    LONG t,*u,*uTemp,*m,*mTemp,*g,*s,*ck,*Wptr,*V;

    cleanArray64s(W,2*(d+1) + 4*(n*d+1) + 2*(dx+1));
    cleanArray64s(Sigmas,n*(d+1));

      for(i=0;i<n;i++){
      for(j=0;j<=d;j++){
         //printf("%d,   ",U[i*(d+1)+j]);
      }
   }

    V = W;
    u = V;       V+= d+1;
    uTemp = V;   V+=  d+1;
    m = V;       V+=  n*d+1;
    mTemp = V;   V+=  n*d+1;
    s = V;       V+=  2*(dx+1);
    ck = V;      V+=  n*d+1;
    g = V;       V+= n*d+1;

    mPos = 0;

    //move c to ck
    dck = 0;
    for (j=0;j<=dx;j++){
       t = c[j];
       ck[j] = t;
       if(t != 0){
          dck = j;
       }

    }

    //printf("c=");polprint64s(ck,dck);printf("\n");

    //Calculate the sigma term and subtract it from C
    for (i=1;i<n;i++){

       cleanArray64s(uTemp,d+1);
       cleanArray64s(mTemp,(n-1)*d+1);
       cleanArray64s(s,2*(dx+1));
       cleanArray64s(g,n*d+1);
       cleanArray64s(V,2*(dx+1));

       //printf("ck:=");polprint64s(ck,dck);printf(":\n");

       //extract M_i,U_i
       point = (d+1)*(i-1);
       for (j=0;j<=d;j++){
          u[j] = U[point + j];
          uTemp[j] = U[point + j];
          if(U[point + j] != 0){
             du = j;
          }
       }

       //printf("ui:= "); polprint64s(uTemp,du); printf(":\n");
       for (j=0;j<(n-i)*d+1;j++){
          m[j] = M[mPos + j];
          mTemp[j] = M[mPos + j];
           if(M[mPos + j] != 0){
              dm = j;
            }
        }
        mPos = mPos + (n-i)*d+1;

       //polgcdext64s(mTemp,uTemp,dm,du,g,s,0,&dg,&ds,&dt,&(W[wTemp]),p);

       //printf("S:");polprint64s(&S[(i-1)*(dx+1)],sDeg[i-1]);printf("\n\n");

       ds = polmul64s(ck,&S[(i-1)*(dx+1)],s,dck,sDeg[i-1],p);
       *numMuls = *numMuls + (dck+1)*(sDeg[i-1]+1);

       //printf("CS:");polprint64s(s,ds);printf("\n\n");

       drem = poldiv64s(s,u,ds,du,p);
       *numMuls = *numMuls + (du)*(ds-du+1);

       //printf("CS div F0:");polprint64s(s,drem);printf("\n\n");

       ds = ds - drem - d;

       for(j=0;j<=drem;j++){
          t = s[j];
          Sigmas[(i-1)*(d)+j] = t;
       }
       for(j=0;j<=ds;j++){
          t = s[drem+j+d];
          s[j] = t;
          s[drem+j+d] = 0;
       }

       //calculate c - r*M_i
       dm = polmul64s(m,&Sigmas[(i-1)*(d)],m,dm,drem,p);
       dm = polsub64s(ck, m, m, dck, dm, p);

       //c is reduced to 0, so the remaining sigmas must be 0.
       if(dm==-1){return;}

       //take quotient of c - r*M_i and U_i;
       drem = poldiv64s(m,u,dm,du,p);
       for(j=0;j<=dm-du;j++){
          t = m[du+j];
          ck[j] = t;
          if(t != 0){
             dck = j;
          }
       }
    }

    //Save final sigmas
    for (j=0;j<=dck;j++){ Sigmas[(n-1)*d+j] = ck[j]; }



    return;
}

//getCoeff2(p1, p2, n, dz , k, k, accDegs,maxDeg, p);
LONG getCoeff2(LONG * Polys, LONG * Coeffs, int n, int dz, int d, LONG p, int start,int *Degs, int *accDeg, LONG *numMuls){

   //local variables
   int i,j,k,t,s,r,MIN,MAX;
   LONG *p1,*p2,*p3;
   ULONG z[2];

   z[0] = 0; z[1] = 0;

   //****************************
   //n=2 case
   //****************************
   if(p < 1LL << 50){
       if(n==2){
          MIN = max32s(0,d-Degs[0]); MAX = min32s(d,Degs[1]);
          k = MIN;
          p1 = &Polys[Degs[0]+1];
          while(k<=MAX-3){
             zfma(z,Polys[k],p1[d-k]);k++;
             zfma(z,Polys[k],p1[d-k]);k++;
             zfma(z,Polys[k],p1[d-k]);k++;
          }
          while(k<=MAX){
             zfma(z,Polys[k],p1[d-k]);k++;
          }
          zmod(z,p);
          return z[0];
        }


       //****************************
       //n>2 case
       //****************************


      //set degree totals
      p1 = &Polys[0];
      p2 = &Polys[Degs[0]+1];
      p3 = &Coeffs[0];

      //initial multiplications
      if(d <= accDeg[1]){
          MIN = max32s(0,d-Degs[1]); MAX = min32s(d,Degs[0]);
          k = MIN;
          while(k<=MAX-3){
             zfma(z,p1[k],p2[d-k]);k++;
             zfma(z,p1[k],p2[d-k]);k++;
             zfma(z,p1[k],p2[d-k]);k++;
          }
          while(k<=MAX){
             zfma(z,p1[k],p2[d-k]);k++;
          }
          zmod(z,p);
          p3[d] = z[0];
      }
       p1 = p2+(Degs[1]+1);

      for (j=1; j<=n-3; j++){

          //make sure i need to do current calculations
          if(d <= accDeg[j+1]){

             //continued polynomial multiplication up to degree d
            MIN = max32s(0,d-accDeg[j]); MAX = min32s(d,Degs[j+1]);
            k = MIN;
            z[0] = 0; z[1] = 0;
            while(k<=MAX-3){
               zfma(z,p1[k],p3[d-k]); k++;
               zfma(z,p1[k],p3[d-k]); k++;
               zfma(z,p1[k],p3[d-k]); k++;
            }
            while(k<=MAX){
               zfma(z,p1[k],p3[d-k]); k++;
            }
            zmod(z,p);
            p3[dz+1+d] = z[0];
         }

         //update pointers
         p1 += Degs[j+1]+1;
         p3 += dz+1;
      }

       //Calculate final coefficient
       MIN = max32s(0,d-accDeg[n-2]);
       MAX = min32s(d,Degs[n-1]);
       z[0] = 0; z[1] = 0;
       j = MIN;
       while(j<=MAX-3){
          zfma(z,p1[j],p3[d-j]);j++;
          zfma(z,p1[j],p3[d-j]);j++;
          zfma(z,p1[j],p3[d-j]);j++;
       }
       while(j<=MAX){
          zfma(z,p1[j],p3[d-j]);j++;
       }
   }
   else{

       if(n==2){
              MIN = max32s(0,d-Degs[0]); MAX = min32s(d,Degs[0]);
              k = MIN;
              p1 = &Polys[dz+1];
              while(k<=MAX-2){
                 zfma(z,Polys[k],p1[d-k]);k++;//if(z[1]>=p){z[1] -= p;}
                 zfma(z,Polys[k],p1[d-k]);k++;if(z[1]>=p){z[1] -= p;}
                 *numMuls = *numMuls + 2;

              }
              while(k<=MAX){
                 zfma(z,Polys[k],p1[d-k]);k++;if(z[1]>=p){z[1] -= p;}
                 *numMuls = *numMuls + 1;
              }
              zmod(z,p);
              return z[0];
            }


           //****************************
           //n>2 case
           //****************************


          //set degree totals
          p1 = &Polys[(n-2)*(dz+1)];
          p2 = &Polys[(n-1)*(dz+1)];
          p3 = &Coeffs[(n-2)*(dz+1)];

          //initial multiplications
          if(d <= accDeg[n-2]){
              MIN = max32s(0,d-Degs[n-1]); MAX = min32s(d,Degs[n-2]);
              k = MIN;
              while(k<=MAX-2){
                 zfma(z,p1[k],p2[d-k]);k++;//if(z[1]>=p){z[1] -= p;}
                 zfma(z,p1[k],p2[d-k]);k++;if(z[1]>=p){z[1] -= p;}
                 *numMuls = *numMuls + 2;
              }
              while(k<=MAX){
                 zfma(z,p1[k],p2[d-k]);k++;if(z[1]>=p){z[1] -= p;}
                 *numMuls = *numMuls + 1;
              }
              zmod(z,p);
              p3[d-(dz+1)] = z[0];
          }
           p1 = p1-(dz+1);
           p3 = p3-(dz+1);

          for (j=n-3; j>=1; j--){

              //make sure i need to do current calculations
              if(d <= accDeg[j]){

                 //continued polynomial multiplication up to degree d
                MIN = max32s(0,d-accDeg[j+1]); MAX = min32s(d,Degs[j]);
                k = MIN;
                z[0] = 0; z[1] = 0;
                while(k<=MAX-2){
                   zfma(z,p1[k],p3[d-k]); k++;//if(z[1]>=p){z[1] -= p;}
                   zfma(z,p1[k],p3[d-k]); k++;if(z[1]>=p){z[1] -= p;}
                   *numMuls = *numMuls + 2;
                }
                while(k<=MAX){
                   zfma(z,p1[k],p3[d-k]); k++;if(z[1]>=p){z[1] -= p;}
                   *numMuls = *numMuls + 1;
                }
                zmod(z,p);
                p3[d-(dz+1)] = z[0];
             }
             p1 = p1-(dz+1);
             p3 = p3-(dz+1);
          }

           //Calculate final coefficient
           MIN = max32s(0,d-accDeg[1]);
           MAX = min32s(d,Degs[1]);
           z[0] = 0; z[1] = 0;
           j = MIN;
           while(j<=MAX-2){
              zfma(z,Polys[j],Coeffs[d-j]);j++;//if(z[1]>=p){z[1] -= p;}
              zfma(z,Polys[j],Coeffs[d-j]);j++;if(z[1]>=p){z[1] -= p;}
              *numMuls = *numMuls + 2;
           }
           while(j<=MAX){
              zfma(z,Polys[j],Coeffs[d-j]);j++;if(z[1]>=p){z[1] -= p;}
              *numMuls = *numMuls + 1;
           }

   }

   zmod(z,p);
   //printf("final solution: %d\n",z[0]);
   return z[0];
}

LONG getCoeff3(LONG * Polys, LONG * Coeffs, LONG * coeffDegs, int n, int dz, int d, LONG p, LONG * shift1, LONG * shift2, LONG * dshift, int *Degs, LONG *W, LONG iter, clock_t *s1, clock_t *s2){

     //denote local variables
     LONG i,j,k,m,t,z,y,s;
     LONG *temp, *tempDegs, *p1, *p2, *p3, *p4, *p5;
     LONG MIN,MAX;
     ULONG Z[2];
     clock_t T1,T2;

    //update calculation
    p1 = Coeffs;
    p2 = Polys;
    for(j=0; j<(n/2); j++){
        if(d <= coeffDegs[j]){
           //printf("got here and j=%d\n",j);
           MIN = max32s(0,d-Degs[2*j+1]);  MAX = min32s(d,Degs[2*j]);
           Z[0] = 0; Z[1] = 0;
           p3 = &p2[Degs[2*j]+1];
           i=MIN;
           while(i <= MAX-3){
              zfma(Z,p2[i],p3[d-i]);i++;
              zfma(Z,p2[i],p3[d-i]);i++;
              zfma(Z,p2[i],p3[d-i]);i++;
           }
           while(i <= MAX){
            zfma(Z,p2[i],p3[d-i]);i++;
           }
           zmod(Z,p);
           p1[shift1[j]+d] = Z[0];
       }
       p2 += coeffDegs[j] + 2;
     }
     if(n%2==1){p1[shift1[n/2]+k] = p2[k];}


     //shift + update coeffcalc main loop
     p2 = p1;
     p1 += dz+n;
     p3 = shift1;
     p5 = coeffDegs;
     t = ceil2(n);

     //update calculation

     //T1 = clock();
     while(t > 1){
        for(j=0; j<(t/2); j++){
          if(d <= p5[n+j]){
               MIN = max32s(0,d-p5[2*j+1]); MAX = min32s(d,p5[2*j]);
               Z[0] = 0; Z[1] = 0;
               p4 = &p2[p5[2*j]+1];
               i = MIN;
               while(i <= MAX-3){
                zfma(Z,p2[i],p4[d-i]); i++;
                zfma(Z,p2[i],p4[d-i]); i++;
                zfma(Z,p2[i],p4[d-i]); i++;
               }
               while(i <= MAX){
                zfma(Z,p2[i],p4[d-i]); i++;
               }
               zmod(Z,p);
               p1[p3[n+j]+d] = Z[0];
          }
          p2 += p5[n+j] + 2;
       }if(t%2==1){p1[p3[n+(t/2)]+d] = p2[d];  }
       p2 = p1;
       p1 += dz+n;
       p3 += n;
       p5+=n;
       t=ceil2(t);
    }

   return p2[d];

}

void createProblem(LONG *A, LONG *F0, int n, int dx, int dz, int dxA, int dzA, LONG alpha, LONG p);

void LagInterpSetup(LONG *A,LONG *xPoints,int dx, LONG *W, LONG p){

    //local variables
    LONG t,a,b,*p1,*p2;
    int i,j,k,m,dm,point1,point2,size;
    ULONG z[2];

    size = (dx+1)/2+1;

    cleanArray64s(W,2*dx+6);
    p1 = &W[dx+2];
    p2 = &W[dx+4];

    //calculate master polynomial
    W[1] = 1;
    p1[1] = 1;
    for(i=1;i<=dx;i++){
       p1[0] = p-xPoints[i];
       polmul64s(W,p1,W,i,1,p);
    }

    //calculate Lagrange polynomials
    for(i=0;i<=dx;i++){
       //copy the master polynomial to spare space
       for(j=0;j<=dx+1;j++){p2[j] = W[j]; }
       p1[0] = p-xPoints[i];
       //perform division
       dm = poldiv64s(p2,p1,dx+1,1,p);
       //move polynomial to A
       for(j=0;j<=dx;j++){ A[i*(dx+1)+j] = p2[j+1]; }

    }

    //calculate inverse and multiply every term by it
    p1 = A;
    for(i=0;i<=dx;i++){
        //calculate denom
        a = 1;
        b = xPoints[i];
        for (j=0;j<=dx;j++){
            if (i!=j){
                a = mul64s(a,sub64s(b,xPoints[j],p),p);
            }
        }
        a = modinv64s(a,p);
        //multiply every term by a
        for(j=0;j<=dx;j++){
            if(p1[j]!=0){
                p1[j] = mul64s(p1[j],a,p);
            }
        }
        p1 = p1 + dx+1;
    }

    //transpose matrix (efficiently)
    cleanArray64s(W,2*(size)*(size-1) + 1);


    if(dx%2 ==0){

        //single entry,
        W[0] = A[0];

        //even case
        p1 = &W[1];
        //firstRow
        for(i=2,k=0;i<=dx+1;i+=2,k++){ p1[k*size] = A[i]; }
        //each column (even elements)
        for(i=2,k=1;i<=dx+1;i+=2,k++){
            for(j=2,m=0;j<=dx+1;j+=2,m++){
                p1[m*(size)+k] = A[(i-1)*(dx+1) + j];
            }
        }

        //odd
        p1 = p1 + size*(size-1);
        //firstRow
        for(i=1,k=0;i<=dx+1;i+=2,k++){ p1[k*size] = A[i]; }
        //each column (even elements)
        for(i=2,k=1;i<=dx+1;i+=2,k++){
            for(j=1,m=0;j<=dx+1;j+=2,m++){
                p1[m*(size)+k] = A[(i-1)*(dx+1) + j];
            }
        }
         //move elements back to A
         for(i=0;i<2*size*(size-1)+1;i++){ A[i] = W[i]; }
    }
    else{
       //transpose matrix A
       for(i=0;i<=dx;i++){
          for(j=i+1;j<=dx;j++){
              t = A[i*(dx+1)+j]; A[i*(dx+1)+j] = A[j*(dx+1)+i]; A[j*(dx+1)+i] = t;
          }
       }
    }

    return;

}

//perform Lagrangian interpolation. The Lagrangian polynomials are stored in LagPolys
//and the interpolation points are stored in yPoints.
//The resulting polynomial is stored in Delta
//The size of W is dx+1
void LagInterpEval(LONG *LagPolys,LONG *yPoints,LONG *Delta, int dx, LONG *W, LONG p, LONG *nummuls){
//Lagpolys - holds the Lagrangian polynomials used for lagrange interpolation - matrix has size (dx+1)^2
//yPoints - array of length dx+1 - contains the y points for Lagrange interpolation (the x points are implied to be 0,1,-1,2,-2,...
//Delta - the interpolated polynomial of degree at most dx
//dx - degree of the main polynomial A in x
//W - additional space needed - size = dx+1
//p - prime

    //local variables
    int i,j,k,size,point;
    LONG *p1,*V;
    ULONG z[2];

    if(dx%2==0){

        size = (dx+1)/2 + 1;
        V = W;
        cleanArray64s(V,size);
        V[0] = yPoints[0];

        //solve first solution (alpha = 0)
        Delta[0] = mul64s(LagPolys[0],yPoints[0],p);
        *nummuls += 1;

        //*******************************************************************
        //solve Lagrange interpolation (even case)
        //*******************************************************************

        //build vector
        for(i=1,j=1;i<size-1;i++,j+=2){
            V[i] = add64s(yPoints[j],yPoints[j+1],p);
        }
        if(dx%2==1){V[size-1] = yPoints[dx];}
        else{V[size-1] = add64s(yPoints[dx-1],yPoints[dx],p);}

        //calculate Delta (even Case)
        p1 = &LagPolys[1];
        if(p < 1LL << 50){
            for(i=0,j=2;i<size-1;i++,j+=2){
               Delta[j] = DotProduct(p1,V,size-1,p);
               p1 = p1+size;
            }
        }
        else{
            for(i=0,j=2;i<size-1;i++,j+=2){
               Delta[j] = DotProductP(p1,V,size-1,p);
               p1 = p1+size;
            }
        }

        //*******************************************************************
        //solve Lagrange interpolation (odd case)
        //*******************************************************************

        //build vector
        for(i=1,j=1;i<size-1;i++,j+=2){
            V[i] = sub64s(yPoints[j],yPoints[j+1],p);
        }
        if(dx%2==1){V[size-1] = yPoints[dx];}
        else{V[size-1] = sub64s(yPoints[dx-1],yPoints[dx],p);}

        //calculate Delta (odd Case)
        if(p < 1LL << 50){
            for(i=0,j=1;i<size-1;i++,j+=2){
               Delta[j] = DotProduct(p1,V,size-1,p);
               p1 = p1+size;
            }
        }
        else{
            for(i=0,j=1;i<size-1;i++,j+=2){
               Delta[j] = DotProductP(p1,V,size-1,p);
               p1 = p1+size;
            }
        }
    }


    //when dx is odd
    else{
        Delta[0] = yPoints[0];
        for(i=1;i<=dx;i++){
           Delta[i] = DotProduct(&LagPolys[i*(dx+1)],yPoints,dx,p);
        }
    }

    return;
}


//Evaluates polynomials at x = 0,+-1,+-2,+-3,...
void FastEval(LONG *Evals, LONG *F, LONG *EvalPoints, int iter, int n, int *degsX, int *degsY, int *degsYprev, int du, int dx, int dz, LONG *W, LONG p){

   int i,j,k,m,ndz,dxdz;
   LONG *p1,*p2,*p3,*p4,e,o,t;
   LONG *shift1, *shift2;
   ULONG z1[2],z2[2];

   ndz = n*(dz+1);
   dxdz = (dx+1)*(dz+1);

   //allocate space for the arrays
   shift1 = W; W += n;
   shift2 = W; W += n;

   //shift elements over to create space.
   shift1[0] = 0;
   shift2[0] = 0;
   for(i=1;i<n;i++){
        shift1[i] = shift1[i-1] + degsY[i-1] + 1;
        shift2[i] = shift2[i-1] + degsYprev[i-1] + 1;
   }

   //perform shift
   for(i=0;i<=dx;i++){
      p1 = &Evals[i*(dz+n)];

      /*printf("before Eval shift\n");
      for(j=0;j<dz+n;j++){
         printf("yyy = %lld ",p1[j]);
      }
      printf("\n");*/

      for(j=n-1;j>=1;j--){
          if(shift1[j] != shift2[j]){
              for(k=degsYprev[j]; k>=0; k--){
                 p1[shift1[j]+k] = p1[shift2[j]+k];
                 p1[shift2[j]+k] = 0;
              }
          }
      }
     /* printf("after Eval shift\n");
      for(j=0;j<dz+n;j++){
         printf("yyy = %lld ",p1[j]);
      }
      printf("\n");*/

   }


   //deal with easy case (x=0, store just constants)
   p1 = &F[iter*(dx+1)];
   for(i=0;i<n;i++){
      if(degsY[i] == iter){Evals[shift1[i] + degsY[i]] = *p1;}
      p1 = p1 + dxdz;
   }

   //preform evaluation if using a smaller prime
   if(p < 1LL << 50){
       p2 = EvalPoints;
       for(i=1;i<=dx/2;i++){
          //pointers
          //p1 = polynomials
          //p2 = x = 0,+-1,...
          //p3 = stores even + odd;
          //p4 = stores even - odd;
          p1 = F + iter*(dx+1);
          p3 = &Evals[((i-1)*2+1)*(n+dz)];
          p4 = &Evals[((i-1)*2+2)*(n+dz)];
          for(j=0;j<n;j++){

             if(degsY[j] == iter){

                 z1[0] = p1[0]; z1[1] = 0;
                 z2[0] = p1[1]; z2[1] = 0;

                 for(m=1,k=2;k<=degsX[j];k+=2,m++){
                    zfma(z1,p1[k],p2[m]);
                    zfma(z2,p1[k+1],p2[m]);
                 }
                 zmod(z1,p); e = z1[0];
                 zmod(z2,p); o = mul64s(z2[0],i,p);

                 p3[ shift1[j] + degsY[j] ] = add64s(e,o,p);
                 p4[ shift1[j] + degsY[j] ] = sub64s(e,o,p);
             }
             p1 = p1 + dxdz;


          }
          p2 = p2 + du+1;

          /*printf("after Eval\n");
          for(j=0;j<dz+n;j++){
          printf("yyy = %lld ",p1[j]);
          }
          printf("\n");*/
       }

       //odd case (if dx%2==1)
       if(dx%2==1){
           p1 = F + iter*(dx+1);
           p3 = &Evals[((dx/2)*2+1)*(n+dz)];
           for(j=0;j<n;j++){
              if(degsY[j] == iter){
                 t=0;e=1;
                 for(k=0;k<=degsX[j];k++){
                    t = add64s(t,mul64s(p1[k],e,p),p);
                    e = mul64s(e,ceil2(dx),p);
                 }

                 p3[ shift1[j] + degsY[j] ] = t;

             }
             p1 = p1 + dxdz;
           }
       }

   }
   else{
      p2 = EvalPoints;
       for(i=1;i<=dx/2;i++){
          //pointers
          //p1 = polynomials
          //p2 = x = 0,+-1,...
          //p3 = stores even + odd;
          //p4 = stores even - odd;
          p1 = F + iter*(dx+1);
          p3 = &Evals[((i-1)*2+1)*(n+dz)];
          p4 = &Evals[((i-1)*2+2)*(n+dz)];
          for(j=0;j<n;j++){

             z1[0] = p1[0]; z1[1] = 0;
             z2[0] = p1[1]; z2[1] = 0;

             for(m=1,k=2;k<=degsX[j];k+=2,m++){
                zfma(z1,p1[k],p2[m]); if(z1[1]>=p){z1[1] -=p;}
                zfma(z2,p1[k+1],p2[m]); if(z2[1]>=p){z2[1] -=p;}
             }
             zmod(z1,p); e = z1[0];
             zmod(z2,p); o = mul64s(z2[0],i,p);

             p3[ shift1[j] + degsY[j] ] = add64s(e,o,p);
             p4[ shift1[j] + degsY[j] ] = sub64s(e,o,p);
             p1 = p1 + dxdz;
          }
          p2 = p2 + du+1;
       }
   }
   return;
}

void CoeffCalcShifts(LONG* CoeffCalc, LONG* Evals, LONG* shift1, LONG* shift2, int* maxDeg, LONG* coeffDegs, LONG n, LONG d, LONG dx, LONG dz, LONG p){

   LONG i,j,k,t,X;
   LONG *p1, *p2, *p3, *p4, *p5;

   //calculate shift positions
      shift1[0]=0; shift2[0]=0;
      for(i=1;i<=n/2;i++){
            shift1[i] = shift1[i-1] + maxDeg[2*(i-1)] + maxDeg[2*(i-1)+1] + 1;
            shift2[i] = shift2[i-1] + coeffDegs[i-1] + 1;
            //coeffDegs[i-1] = z;
      }
      t = ceil2(n);
      p1 = shift1 + n;   p2 = shift2 + n;
      p3 = shift1; p4 = coeffDegs + n;
      while(t>1){
         p1[0] = 0; p2[0] = 0;
         for(i=1;i<=t/2;i++){
            p1[i] = p1[i-1] + p3[2*(i-1)+2] - p3[2*(i-1)] - 1;
            p2[i] = p2[i-1] + p4[i-1] + 1;
        }

        p1 += n; p2 += n; p3 += n; p4 += n;
        t = ceil2(t);
     }

     //perform shift
     //perform initial shift and calculation
     for(i=0;i<=dx;i++){

         p1 = &CoeffCalc[i*log2N(n)*(dz+n)];
         p2 = &Evals[i*(dz+n)];

         //shift
         if(n%2==1){ //if there's an odd # of elements
             for(k=maxDeg[n-1];k>=0;k--){
                p1[shift1[n/2]+k] = p1[shift2[n/2] + k];
                if(shift1[n/2] != shift2[n/2]){p1[shift2[n/2] + k] = 0;}
             }
         }
         t = (n/2) - 1;
         for(j=t;j>=1;j--){
             for(k=maxDeg[2*j]+maxDeg[2*j+1];k>=0;k--){
                 p1[shift1[j]+k] = p1[shift2[j]+k];
                 if(shift1[j] != shift2[j]) {p1[shift2[j]+k] = 0;}
             }
         }


         p3 = shift1;
         p4 = shift2;
         p5 = coeffDegs;
         t = ceil2(n);

         while(t>1){

            p1 += dz+n;
            p3 += n; p4 += n;

            //preform standard shift
            if(t%2==1){ //if there's an odd # of elements
               if(p3[t/2] != p4[t/2]){
                 for(k=p5[n+(t/2)-1];k>=0;k--){
                    p1[p3[t/2]+k] = p1[p4[t/2] + k];
                    p1[p4[t/2] + k] = 0;
                 }
              }
          }
            for(j=(t/2)-1;j>=1;j--){
                if(p3[j] != p4[j]){
                    for(k=p5[n+j];k>=0;k--){
                        p1[p3[j]+k] = p1[p4[j]+k];
                        p1[p4[j]+k]=0;
                    }
                }
            }

            t = ceil2(t);
         }



     }
    return;
}

//This preforms hensel lifting on n factors to factor polynomial A
int HenselLiftCubic(LONG *A, int dx, int dz, LONG *F0, int n, int du, LONG *F, LONG alpha, LONG *TempSpace, LONG p ){
/*
A[x][z]- polynomial you wish to factor
dx - degree bound on variable x of A
dz - degree bound on variable z of A
F0 - set of linear factors
n - number of initial factors
du - degree bound of initial factors (could be extended to an array later)
F - array that will store the final polynomials. must be of size n*(dx+1)*(dz+1)
W - working array for calculations. Must be at least size (tbd)
alpha - LONG integer
p - prime (duh)
*/

   //local variables
   int i,j,k,l,maxSizeA,ndz,flag,*maxDeg, *maxDegPrev,*sDeg, *accDegs, *curXDeg, sumDeg, InterpFlag;
   LONG t,t2,z,*M,*E,*p1,*p2,*p3,*p4,*p5,*TEMP1,*Coeffs, *LagInterpPolys;
   LONG *DioSVal, numDioMuls,numCoeffExtractMul,numEvalMuls,numInterpMuls, * coeffDegs, *dshift, *shift1, *shift2, *g, *fac, *defac;
   LONG *CoeffCalc,*Evals,*ck,*evalPoints,*interpPoints,*Delta,*EvalPointsMul, *tempDegs,*W, wTemp;
   LONG temp,temp2,temp3;
   recint P;
   clock_t T1,T2,T3,T4,s1,s2,s3,s4,s5,s6,s7,s8,s9,s10,s11,s12,s13,s14,s15;

   //set time to 0
   s1=0;s2=0;s3=0;s4=0;s5=0;s6=0;s7=0;s8=0;s9=0;s10=0;s11=0;s12=0;s13=0;s14=0;s15=0;

   //initial size calculations
   maxSizeA = (dx+1)*(dz+1);
   ndz = n*(dz+1);
   numEvalMuls = 0;
   numDioMuls = 0;
   numCoeffExtractMul = 0;
   numInterpMuls = 0;
   InterpFlag = 0;

   P = recip1(p);

   //declare arrays
   W = TempSpace;
   maxDeg = arrayI(n); //this is in y
   maxDegPrev = arrayI(n);
   accDegs = arrayI(n); //accumulative sum of degrees in y (for coeffcalc)
   sDeg = arrayI(n);  //Degree of each s variable (for Diophantine Equation)
   curXDeg = arrayI(n);  //Degree of x for each factor (efficiency)
   M = W;                          W += (n-1)*n*du/2 + n - 1; //M polynomials (diophantine)
   E = W;                          W += maxSizeA;             //The error at each iteration
   TEMP1 = W;                      W += maxSizeA;             //Temporary storage
   evalPoints = W;                 W += dx+1;
   interpPoints = W;               W += dx+1;
   LagInterpPolys = W;             W += (dx+1)*(dx+1);
   Delta = W;                      W += dx+1;
   CoeffCalc = W;                  W += (dx+1)*log2N(n)*(dz+n);
   Evals = W;                      W += (dx+1)*(dz+n);
   EvalPointsMul = W;              W += (du+1)*((dx/2)+1);
   DioSVal = W;                    W += (n-1)*(dx+1);
   coeffDegs = W;                  W += log2N(n)*n;
   tempDegs = W;                   W += log2N(n)*n;
   dshift = W;                     W += n;
   shift1 = W;                     W += n*(log2N(n)+1);
   shift2 = W;                     W += n*(log2N(n)+1);
   g = W;                          W += dz+1;
   fac = W;                        W += dz+1;
   defac = W;                      W += dz+1;

   //move contents of F0 to F
   //Calculate exact degree bound of factors
   p1 = F;
   p2 = F0;
   for(i=0;i<n*log2N(n);i++){coeffDegs[i]=-1;}
   cleanArray64s(shift1,log2N(n)*n);
   cleanArray64s(shift2,log2N(n)*n);
   cleanArray64s(tempDegs,log2N(n)*n);

   for(i=0;i<n;i++){
      curXDeg[i] = -1;
      for(j=0;j<=du;j++){
         if(p2[j] != 0){
            curXDeg[i]=j;
            p1[j] = p2[j];
         }
      }
      p1 = p1 + maxSizeA;
      p2 = p2 + du+1;
   }

   //declare maxDeg, the size of each n factors
   for(i=0;i<n;i++){
      maxDeg[i] = 0;
      maxDegPrev[i] = 0;
   }

   //set evalPoints
   evalPoints[0] = 0;
   p1 = evalPoints+1;
   for(i=1;i<=dx/2;i++){
      p1[0] = i;
      p1[1] = p-i;
      p1 +=2;
   }
   if(dx%2 == 1){
      evalPoints[dx] = i;
   }

   //setup shift1 + shift2
   for(i=0;i<=n/2;i++){
      shift1[i]=i;
      shift2[i]=i;
   }
   t = ceil2(n);
   p1 = shift1+n; p2 = shift2+n;
   while(t>1){
    for(i=0;i<=t/2;i++){
       p1[i]=i; p2[i]=i;
    }
    p1 += n; p2 += n;
    t=ceil2(t);
   }

   //Generate Lagrange Polynomials for Interpolation in main loop
   LagInterpSetup(LagInterpPolys,evalPoints,dx,W,p);

   //calculate evaluation point constants
   p1 = EvalPointsMul;
   for(i=1;i<=dx/2;i++){
      t = 1;
      p1[0] = t;
      t2 = i*i;
      for(j=1;j<=du;j++){
         t = mul64s(t,t2,p);
         p1[j] = t;
      }
      p1 = p1 + du+1;
   }

   //Recover multiple of U for future calls to diophantine equation
   j = MultiEEA(F0,du,n,M,W,p);

   //if the EEA failed, return 'FAIL'(add last)
   if(j==-1){
      printf("The input factors are not coprime");
      return -1;
   }

   //Generate the S values from Gcdex
   generateGcdexS(F0,n,du,M,dx,DioSVal,sDeg,W,p);

   // declare initial Error
   for(i=0;i<maxSizeA;i++){
      E[i] = A[i];
   }
   //Change it to its Taylor Representation E = a_0 + a_1(y-alpha) + ...
   T1 = clock();
   if(alpha != 0){altDiv2(E,dx,dz,alpha,W,p);}
   s9 = clock()-T1;

   //printf("E:=");polprintbivar64s(E,dx,dz,"x","(y-5)");printf(":\n");

   //Preform initial polynomial evaluations
   FastEval(Evals, F,EvalPointsMul, 0, n, curXDeg, maxDeg, maxDegPrev, du, dx, dz, W, p);

   //print out evals 2

   /*temp = 0;
   for(i=0;i<=dx;i++){
       printf("x=%d\n",i);
       for(j=0;j<dz+n;j++){
          printf("%d ",Evals[temp+j]);
       }
       printf("\n");
       temp = (i+1)*(dz+n);
   }*/

   //initial CoeffCalc calculations
   p1 = CoeffCalc;
   p2 = Evals;
   for(i=0;i<=dx;i++){
      for(j=0;j<n/2;j++){ p1[j] = mul64s(p2[2*j],p2[2*j+1],p); }
      if(n%2==1){ p1[n/2] = p2[n-1]; }
      p3 = p1;
      t=ceil2(n);

      while(t>1){
        for(j=0;j<t/2;j++){p3[(dz+n) + j] = mul64s(p3[2*j],p3[2*j+1],p);}
        if(t%2==1){ p3[(dz+n) + (t/2)] = p3[t-1]; }
        t = ceil2(t);
        p3 += dz+n;
      }


      p1 += log2N(n)*(dz+n);
      p2 += (n+dz);
   }

    //print out coeffcalc

        /*printf("x=1\n");
        for(j=0;j<log2N(n);j++){
            printf("n=%d\n",j);
            for(k=0;k<dz+n;k++){
                printf("z=%d ",CoeffCalc[log2N(n)*(dz+n) + j*(dz+n) + k]);
            }
            printf("\n");
        }*/


   //main FOR loop
   T3 = clock();
   ck = E;
   for(k=1;k<=dz;k++){

        //printf("Iteration: k = %d\n\n",k);

      //get the error for this loop
      ck += dx+1;


     /*printf("shift1 info\n");
    for(i=0;i<log2N(n);i++){
        for(j=0;j<n;j++){
            printf("%lld ",shift1[i*n + j]);
        }
        printf("\n");
    }
    printf("shift2 info\n");
    for(i=0;i<log2N(n);i++){
        for(j=0;j<n;j++){
            printf("%lld ",shift2[i*n + j]);
        }
        printf("\n");
    }
    printf("coeffDeg info\n");
    for(i=0;i<log2N(n);i++){
        for(j=0;j<n;j++){
            printf("%lld ",tempDegs[i*n + j]);
        }
        printf("\n");
    }*/


         //Preform coefficient extraction
         T1 = clock();
          p1 = Evals;
          p2 = CoeffCalc;



          for(i=0;i<=dx;i++){
             //calculate coefficient for (y-alpha)^k
             interpPoints[i] = getCoeff3(p1,p2,tempDegs,n,dz,k,p,shift1,shift2,dshift,maxDeg,W,i,&s14,&s15);
             p1 += (n+dz);  p2 += log2N(n)*(dz+n);
             //printf("zzz = %lld ",interpPoints[i]);
          }
          //printf("\n");
          T2 = clock();
          s2 = s2 + T2 - T1;


          /*printf("after getCoeff3\n");
          //for(t=0;t<=dx;t++){
          for(i=0;i<log2N(n);i++){
            for(j=0;j<dz+n;j++){
                printf("z = %lld ",CoeffCalc[log2N(n)*(dz+n) + i*(dz+n) + j]);
            }
            printf("\n");
          }
          printf("\n");*/
          //}

          //update CoeffDeg
          for(i=0;i<n/2;i++){ coeffDegs[i] = maxDeg[2*i] + maxDeg[2*i+1]; }
          if(n % 2 == 1){ coeffDegs[n/2]= maxDeg[n-1];}
          t = ceil2(n);
          p1 = coeffDegs;
          while(t>1){
             for(j=0;j<t/2;j++){p1[n + j] = p1[2*j] + p1[2*j+1];}
             if(t%2 == 1){p1[n + (t/2)] = p1[t-1];}
             p1 += n;
             t = ceil2(t);
          }


      //Interpolate
      cleanArray64s(Delta,dx+1);
      LagInterpEval(LagInterpPolys,interpPoints,Delta,dx, W, p, &numInterpMuls);

      //printf("Delta := "); polprint64s(Delta,dx); printf("\n");

      //printf("temp=");polprint64s(ck,dx);printf("\n");

      //subtract delta from e
      temp = polsub64s(ck, Delta, ck, dx, dx, p);

      //printf("temp=");polprint64s(ck,temp);printf("\n");

      //check to make sure ck isn't zero
      flag = 0;
      for (i=0;i<=dx;i++){
        if (ck[i] != 0){
            flag = 1;
            break;
        }
      }

      //add the sum of all the degrees
      sumDeg = 0;
      for (i=0;i<n;i++){
         sumDeg += maxDeg[i];
      }

      //check fail case
      if(flag == 1 && sumDeg == dz){
         printf("Error: c_k != 0 when the sum of degrees is equal to d\n\n");
         return -1;
      }
      //if c_k isn't equal to zero, solve diophantine equation and preform updates/evaluations
      else if(flag == 1){
          //Solve Diophantine Equation
          T1 = clock();

          cleanArray64s(TEMP1,maxSizeA);

          Dionvar(F0, n, du, ck, dx, M, DioSVal,sDeg, W, TEMP1, p, &numDioMuls);
          T2 = clock();
          s4 = s4 + T2 - T1;


          /*for(i=0;i<n;i++){
             polprint64s(&TEMP1[i*(du)],du-1);printf("\n");
          }
          printf("\n\n");*/

          //Update F
          //Update degree of X counter for factors
          p1 = &F[k*(dx+1)];
          p2 = TEMP1;
          InterpFlag = 1;
          for(i=0;i<n;i++){
             curXDeg[i]=-1;
             for(j=0;j<du;j++){
                if(p2[j] != 0){
                   InterpFlag = 0;
                   curXDeg[i]=j;
                   maxDeg[i] = k;
                   p1[j] = p2[j];
                }
             }
             p1 = p1 + maxSizeA;
             p2 = p2 + du;
          }

          for(i=0;i<n;i++){
            //printf("f%d:=",i+1);polprintbivar64s(&(F[i*(dx+1)*(dz+1)]),dx,dz,"x","(y-5)"); printf(":\n");
          }

          //generate accumulted degrees of polynomials (for coefficient extraction
          accDegs[0] = maxDeg[0];
          for(i=1;i<n;i++){
             accDegs[i] = accDegs[i-1] + maxDeg[i];
          }

          /*printf("just checking degrees\n");
          for(i=0;i<n;i++){
            printf("n=%d n=%d\n",maxDeg[i],maxDegPrev[i]);
          }
          printf("\n");*/

          //Calculate evaluation points
          FastEval(Evals, F,EvalPointsMul, k, n, curXDeg, maxDeg, maxDegPrev, du, dx, dz, W, p);


          //print Evals
          /*for(j=0;j<=dx;j++){
              for(i=0;i<dz+n;i++){
                printf("yyy = %lld ",Evals[j*(dz+n)+i]);
              }
          printf("\n");
          }
          printf("\n");*/


          //update degrees
          for(i=0;i<n;i++){maxDegPrev[i] = maxDeg[i];}

          //calculate shifts for getcoeff3
          dshift[0]=0;
          for(i=1;i<n;i++){
            dshift[i] = dshift[i-1] + maxDeg[i-1] + 1;
          }

          /*printf("before shift\n");
          //for(t=0;t<=dx;t++){
          for(i=0;i<log2N(n);i++){
            for(j=0;j<dz+n;j++){
                printf("z = %lld ",CoeffCalc[log2N(n)*(dz+n) + i*(dz+n) + j]);
            }
            printf("\n");
          }
          printf("\n");*/
          //}

          //perform shift
          CoeffCalcShifts(CoeffCalc,Evals,shift1,shift2,maxDeg,coeffDegs,n,k,dx,dz,p);


          /*printf("after shift\n");
          for(i=0;i<log2N(n);i++){
            for(j=0;j<dz+n;j++){
                printf("z = %lld ",CoeffCalc[log2N(n)*(dz+n) + i*(dz+n) + j]);
            }
            printf("\n");
          }
          printf("\n");

          printf("just checking degrees\n");
          for(i=0;i<n;i++){
            printf("n=%d ",maxDeg[i]);
          }
          printf("\n");*/

          //Update CoeffCalc
          cleanArray64s(TEMP1, n);
          for(i=0;i<=dx;i++){
             p1 = &Evals[i*(dz+n)];
             p2 = &CoeffCalc[i*log2N(n)*(dz+n)];
             p5 = tempDegs;
             for(j=0;j<n/2;j++){
                TEMP1[j] = 0;
                if(maxDeg[2*j]==k) {TEMP1[j] = mul64s(p1[k],p1[maxDeg[2*j]+1],p);}
                if(maxDeg[2*j+1]==k) {TEMP1[j] = add64s(TEMP1[j],mul64s(p1[0],p1[maxDeg[2*j]+1+k],p),p);}
                p2[k] = add64s(p2[k],TEMP1[j],p);
                tempDegs[j] = maxDeg[2*j] + maxDeg[2*j+1];
                p1 += maxDeg[2*j] + maxDeg[2*j+1] + 2;
                p2 += maxDeg[2*j] + maxDeg[2*j+1] + 1;
             }
             if(n%2==1){p2[k] = p1[k]; TEMP1[n/2] = p1[k]; tempDegs[n/2]=maxDeg[n-1];}
             p2 = &CoeffCalc[i*log2N(n)*(dz+n)];
             p3 = coeffDegs;
             p4 = shift1;
             t = ceil2(n);
             //if(i==1){for(j=0;j<n;j++){ printf("s=%lld ",TEMP1[j]); }printf("\n");}
             while(t>1){
                for(j=0;j<t/2;j++){
                   if(p5[2*j]>=k){ TEMP1[j] = mul64s(TEMP1[2*j],p2[p4[2*j+1]],p); } else {TEMP1[j]=0;}
                   if(p5[2*j+1]>=k){TEMP1[j] = add64s(TEMP1[j],mul64s(TEMP1[2*j+1],p2[p4[2*j]],p),p);  }
                   p2[(dz+n)+ p4[n+j]+k] = add64s(p2[(dz+n)+ p4[n+j]+k],TEMP1[j],p);
                   p5[n+j] = p5[2*j]+p5[2*j+1];
                }
                if(t%2==1){p2[(dz+n) + p4[n+(t/2)] + k] = add64s(p2[(dz+n) + p4[n+(t/2)] + k],TEMP1[t-1],p); TEMP1[(t/2)] = TEMP1[t-1]; p5[n+(t/2)] = p5[t-1]; }
                p2 += dz+n; p3 += n; p4 += n; p5 += n;
                t = ceil2(t);
                //if(i==1){for(j=0;j<n;j++){ printf("s=%lld ",TEMP1[j]); }printf("\n");}
             }
         }
         /*printf("after Coeff Update\n");
          for(i=0;i<log2N(n);i++){
            for(j=0;j<dz+n;j++){
                printf("z = %lld ",CoeffCalc[log2N(n)*(dz+n) + i*(dz+n) + j]);
            }
            printf("\n");
          }
          printf("\n");*/

      }

   }

   T4 = clock();
   cleanArray64s(TEMP1,dx+1);

   //check if sum of degrees = dzvpn.its.sfu.ca
   if(sumDeg < dz){
      printf("Error: sum of degrees not equal to %d\n",dz);
      return -1;
   }

   //print statements
   //printf("Change of Base time=%8.2fms\n",(s9)/1000.0/1.0);
   //printf("Evaluation time=%8.2fms\n",(s1)/1000.0/1.0);
   //printf("Coefficient Extraction time=%8.2fms\n",(s2)/1000.0/1.0);
   //printf("Interpolation Setup time=%8.2fms\n",(s7)/1000.0/1.0);
   //printf("Interpolation Calc time=%8.2fms\n",(s3)/1000.0/1.0);
   //printf("Interpolation Calc part 1 time=%8.2fms\n",(s14)/1000.0/1.0);
   //printf("Interpolation Calc part 2 time=%8.2fms\n",(s15)/1000.0/1.0);
   //printf("Total Interpolation time=%8.2fms\n",(s3+s7)/1000.0/1.0);
   //printf("Diophantine time=%8.2fms\n",(s4)/1000.0/1.0);
   /*printf("The number of Evaluation Multiplications is %lld\n",numEvalMuls);
   printf("The number of Diophantine Multiplications is %lld\n",numDioMuls);
   printf("The number of Interpolation Multiplications is %lld\n",numInterpMuls);
   printf("The number of Coefficient Extraction Multiplications is %lld\n",numCoeffExtractMul);*/
   //printf("\n");


   free(maxDeg); free(accDegs); free(sDeg); free(curXDeg);

   return 0;

}


// a(x) = a0 + a1 x + a2 x^2 + a3 x^3 + a4 x^4  in Zp[x]
// a(x+b) = a0 + a1 (x+b) + a2 (x+b)^2 + a3 (x+b)^3 + a4 (x+b)^4  in Zp[x]
//        = a0 + (x+b) [ a1 + (x+b) [ a2 + (x+b) [ a3 + (x+b) a4 ] ] ]

LONG TaylorShift2( LONG *a, LONG d, LONG b, LONG p ) {
// a(x+b) using Horner's method (in place)
    LONG i, k;
    if( d==-1 ) return -1;
    if( d==0 ) return 0; 
    for ( k=d-1; k>=0; k-- ) {
        for ( i=k; i<=d-1; i++ ) 
            a[i] = add64s( a[i],mul64s(b,a[i+1],p),p );
    }
    return d;
}

void ExpandF( LONG *F, int r, LONG dx, LONG dy, LONG alpha, LONG p) {
	// F is the answer of BHL (HenselLiftCubic), where F[rr][(dx+1)*j+i] = coeff( coeff(f[rr],y-5,j), x, i)
	// The input alpha is p - alpha[j] in CMSHL Maple code
	// Return expanded f onto F
	LONG *Temp=array(dy+1), i, j, k, d; int rr; 
	for ( rr=1; rr<=r; rr++ ) {
	    for ( i=0; i<dx+1; i++ ) {
		    for ( j=0; j<dy+1; j++ ) Temp[j] = F[(rr-1)*(dx+1)*(dy+1) + (dx+1)*j+i]; 
			k = dy; while ( k>=0 ) { if ( Temp[k] == 0 ) k--; else break; } 
			d = k;
		    TaylorShift2( Temp, d, alpha, p ); 
	        for ( j=0; j<dy+1; j++ ) F[(rr-1)*(dx+1)*(dy+1) + (dx+1)*j+i] = Temp[j];
	    }
	}
	free(Temp);
	return;
}

//main
/* int main () {

    clock_t ST,T1,T2,T3,T4;
    LONG p;

    int dx,dz,n,du,i,j,k,s,point1,point2,db1,db2,size;
    LONG *A,*B,*C,*D,*F0,*F,*W,*disp,alpha,t,*Z1,*Z2,x,y,tempArraySize;
    recint P;
    ULONG z[2];

   seed = 1;
   mult = 6364136223846793003ll;
   p = pow(2,31)-1;
   //p = 101;du  = du*2;
   //p = 1009;
   //p = 2179;zfma(Z,p2[i],p4[d-i]); i++;
   //p = 1125899906842597;
   //p = pow(2,51) + 21;
   //p = 1152921504606846883;
   //du = 3;
   P = recip1(p);

    //for(s=1;s<13;s++){

        n = 3;
        //highest degree of factors
        dx = 7;
        dz = 5;
        du = 3;
        alpha = 5;
        A = array((dx+1)*(dz+1));
        disp = array((dx+1)*(dz+1));
        F0 = array((du+1)*n);
        F = array((dx+1)*(dz+1)*n);
        tempArraySize = n + n*(n-1)*dx + 2*(dx+1)*(dz+1) + 9*(dx+1) + 4*(dx+1)*(dx+1) + 2*n*(dx+1)*(dz+1) + 5*n*(dx+1) + 32*(dz+1);
        //tempArraySize = (dx/n)*(4*n+1) + 4*dx + 10 + 50*(dx+1)*(dz+1);
        W = array(tempArraySize);

        //printf("There are %d factors, each of which has degree(x) = %d and degree(z) = %d\n",n,du,du);
        printf("The total degree of A is degree(x) = %d and degree(z) = %d\n",dx,dz);

        cleanArray64s(A,(dx+1)*(dz+1));
        cleanArray64s(F,(dx+1)*(dz+1)*n);
        cleanArray64s(W,tempArraySize);
        cleanArray64s(F0,(du+1)*n);
        cleanArray64s(disp,(dx+1)*(dz+1));

        A[0] = 0; A[1] = 70; A[2] = 2147483427; A[3] = 53; A[4] = 2147483449; A[5] = 2147483631; A[6] = 22; A[7] = 1;
	A[(dx+1)] = 2147483213; A[(dx+1)+1] = 1984; A[(dx+1)+2] = 2147481625; A[(dx+1)+3] = 3602;
        A[(dx+1)+4] = 2147482468; A[(dx+1)+5] = 1706; A[(dx+1)+6] = 167;

        A[2*(dx+1)] = 2147479803; A[2*(dx+1)+1] = 8176; A[2*(dx+1)+2] = 2147473501; A[2*(dx+1)+3] = 14827;
        A[2*(dx+1)+4] = 2147478446; A[2*(dx+1)+5] = 6860;

        A[3*(dx+1)] = 2147476765; A[3*(dx+1)+1] = 8682; A[3*(dx+1)+2] = 2147466499; A[3*(dx+1)+3] = 9280;
        A[3*(dx+1)+4] = 2147475339;

	A[4*(dx+1)] = 2147480175; A[4*(dx+1)+1] = 4884; A[4*(dx+1)+2] = 2147475339; A[4*(dx+1)+3] = 2464;
        A[5*(dx+1)+1] = 2464;
      //  A[5*(dx+1)] = 2147479671;

	F0[0] = 2147481965; F0[1] = 457; F0[2] = 1; F0[3] = 0;
	F0[4] = 310; F0[5] = 2147482537; F0[6] = 400; F0[7] = 1;
        F0[8] = 6; F0[10] = 1;





        printf("A := ");polprintbivar64s(A,dx,dz,"x","y");
        printf("\n\n");
        for(i=0;i<n;i++){
           printf("f%d :=",i+1);polprint64s(&F0[i*(du+1)],du); printf(":");
           printf("\n");
        }
        //printf("%d\n\n",dx/n);
        T1 = 0;
        size = 1;
        for(i=0;i<size;i++){
            cleanArray64s(F,(dx+1)*(dz+1)*n);
            cleanArray64s(W,tempArraySize);
            T3 = clock();
            HenselLiftCubic(A, dx, dz, F0, n, du, F, alpha, W, p);
            T4 = clock();
            T1 = T1 + (T4-T3);
        }
        printf("Time for Cubic Algorithm = %8.2fms\n",T1/1000.0/1.0/size);
        printf("\n");

        printf("F := \n");
        for(i=0;i<n;i++){
            polprintbivar64s(&F[i*(dx+1)*(dz+1)],dx,dz,"x","(y-5)");
            printf("\n\n");
        }

        //polprintbivar64s(F,dx,n*(dz+1),"x","y");
        printf("Comparison: \n");
        ComparePolyFactorsBivar(A,F,dx,dz,n,alpha,p);
        printf("\n");

        free(A);
        free(disp);
        free(F0);
        free(F);
        free(W);

    //}
    return 0;

} */