//Author-Ayoola Jinadu
#include <stdio.h>
#include <stdlib.h>

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

LONG addmod(LONG a,LONG b,LONG p)  // a+b mod p
{   // assume 0 <= a < p, 0 <= b < p
    LONG s;
    s = (a-p)+b; // -p <= s <= p-2
    if( s<0 ) return s+p; else return s;
}

LONG submod1(LONG a, LONG b, LONG p )
{  LONG t;
   t = a-b;
   if( t<0 ) return t+p; else return t;
}

LONG submod(LONG a, LONG b, LONG p) { LONG t; t = a-b; t += (t>>63) & p; return t; }

LONG min(LONG a, LONG b) { if( a<b ) return a; else return b; }
LONG max(LONG a, LONG b) { if( a>b ) return a; else return b; }
LONG neg(LONG a, LONG p) { return (a==0) ? 0 : p-a; }



LONG mulmod(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;
}



/* c^(-1) mod p assuming 0 < c < p < 2^63 */
/* if gcd(c,p) <> 1 then it returns 0 */
LONG invmod( LONG c, LONG p )
{   LONG d,r,q,r1,c1,d1;
    d = p;
    c1 = 1;
    d1 = 0;
    while( d != 0 ) {
        q = c / d;
        r = c - q*d; r1 = c1 - q*d1;
        c = d; c1 = d1;
        d = r; d1 = r1;
    }
    if( c!=1 ) return( 0 );
    if( c1 < 0 ) c1 += p;
    return( c1 );
}

LONG powmod(LONG a, LONG n, LONG p )
{  // a^n mod p where 0 <= a < p
   LONG t;
   if( a==0 ) return 0;
   if( n==0 ) return 1;
   //  a^9 mod p  =   a(a^4)^2  all mod p
   t = powmod(a,n/2,p);
   t = mulmod(t,t,p);
   if( n&1 ) t = mulmod(t,a,p);
   return t;
}

ULONG seed;
ULONG mult;

LONG randmod(LONG p) {
    // output x in [0,p-1]
    LONG x,y;
    extern ULONG seed, mult;
    seed = mult*seed;
    x = seed >> 32;
    seed = mult*seed;
    y = seed >> 32;
    x = (x<<31) | y;
    x = x % p;
    return(x);
}


/***********************************************************************/
/*                Array routines                                       */
/***********************************************************************/
LONG *array( LONG n )
{  LONG *A;
if( n>1000 ) printf("n=%lld\n",n);
   A = (LONG *) malloc( 8*n );
   if( A==0 ) { printf("out of memory\n"); exit(1); }
   return A;
}

void arrayprint( LONG *A, int n )
{   int i;
    printf("[");
    for( i=0; i<n; i++ ) { printf("%lld",A[i]); if( i<n-1 ) printf(", "); }
    printf("];\n");
    return;
}

/***********************************************************************/
/*                Polynomial routines                                  */
/***********************************************************************/

void printpol(LONG *f, LONG d) {LONG i;
  if( d==-1 ) {printf("0"); return;}
  if( f==0 ) {printf("0"); return;}
    for( i=d; i>0; i-- ){ if( f[i]!=0 ) printf("%lld*y^%lld+", f[i], i );}
    printf("%lld", f[0]);
  return;
}

int polequal(LONG *a, LONG da, LONG *b, LONG db) {
   LONG i;
   if( da<0 && db<0 ) return 1;
   if( da!=db ) return 0;
   for( i=0; i<=da; i++ ) if( a[i]!=b[i] ) return 0;
   return 1;
}

//if the 2 polynomials are the same return 1


LONG poladd( LONG *a, LONG da, LONG *b, LONG db, LONG *c,  LONG p )
{  // c = a + b in Fp[x]  and c is space for the answer assumed big enough
   LONG i,dc;
   if( da>db ) return poladd(b,db,a,da,c,p);
   for( i=0; i<=da; i++ ) c[i] = addmod(a[i],b[i],p);
   for( i=da+1; i<=db; i++ ) c[i] = b[i];
   // [1,2,3] + [4,5,8] mod 11 = [5,7,0]  return 1
    for( dc=db; dc>=0 && c[dc]==0; dc-- );
    return dc; //Here its returning a LONG and c is pointing to the address of that long
}

LONG polsub(LONG *a, LONG da,LONG *b, LONG db, LONG *c,  LONG p) {
// c = a-b mod p
   LONG i,m;
   m = min(da,db);
   for( i=0; i<=m; i++ ) c[i] = submod(a[i],b[i],p);
   if( da==db ) { while( da>=0 && c[da]==0 ) da--; return da; }
   if( da>db ) { for ( i=db+1; i<=da; i++ ) c[i] = a[i]; return da; }
   else { for( i=da+1; i<=db; i++ ) c[i] = neg(b[i],p); return db; }
}


LONG  polmul( LONG *a, LONG m, LONG *b, LONG n, LONG *c, LONG p )
{  // c = a * b in Fp[x]  and c is space for the answer assumed big enough
   // (a0+a1*x+...+am*x^m) x (b0+b1*x+...+bn*x^n)
    LONG i,j,dc;
   if( m==-1 || n==-1 ) return -1;
   for( i=0; i<=m+n; i++ ) c[i] = 0;
   for( i=0; i<=m; i++ )
      for( j=0; j<=n; j++ ) c[i+j] = addmod( c[i+j], mulmod(a[i],b[j],p), p );
   for( dc=m+n; dc>=0 && c[dc]==0; dc-- )
       ;
    return dc;
    
}



LONG * treemul1( LONG *S, LONG n, LONG p ) {
   LONG m, *Q, *Q1, *Q2;
   Q = array(n+1);
   if( n==0 ) { Q[0] = 1; return Q; }
   if( n==1 ) { Q[0] = submod(0,S[0],p); Q[1] = 1; return Q; }
   m = n/2;
   Q1 = treemul1( S, m, p );
   Q2 = treemul1( S+m, n-m, p );
   polmul( Q1, m, Q2, n-m, Q, p );
   free(Q1);free(Q2);
   return Q;
}

LONG horner(LONG *a, LONG d, LONG x, LONG p) {
   LONG r,i;
   if( d==-1  ) return 0;
   r = a[d];
   if (d==0)   return r;
   for( i=d-1; i>=0; i-- )
      r = addmod( mulmod(r,x,p), a[i], p );
   return r; // We do d mults and d adds
}

LONG poldiv( LONG *a, LONG da, LONG *b, LONG db, LONG *q, LONG p )
{   //  Compute a div b so that a = b * q + r  with quotient q and remainder r
    //  Output the remainder in a of degree dr and quotient in q with deg(q) = n-m
    //  Allow q to be 0 meaning we only compute the remainder -- for gcd
    LONG i,dr,inv,t;
    if( db==-1 ) { printf("division by 0\n"); exit(1); }
    if( da<db ) return da;
    inv = invmod(b[db],p);
    for( dr=da; dr>=db; dr-- ) {
        if( a[dr]!=0 ) { // a = a - a[dr]/b[db] x^(dr-db) b
            t = mulmod(a[dr],inv,p);
           if( q!=0 )  q[dr-db] = t;
            for( i=db-1; i>=0; i-- )
                a[dr-(db-i)] = submod(a[dr-(db-i)],mulmod(t,b[i],p),p);
            a[dr] = 0;
        }
        else if( q!=0 ) q[dr-db] = 0;
    }
    while( dr>=0 && a[dr]==0 ) dr--;
    return dr;
}



LONG polscamul( LONG *a, LONG d, LONG z, LONG p ) {
    // a(x) = z a(x) mod p
    LONG i;
    for( i=0; i<=d; i++ ) a[i] = mulmod(a[i],z,p);
    while( d>=0 && a[d]==0 ) d--;
    return d;
}

LONG polmonic( LONG *a, LONG d, LONG p ) {
    if( d<0 ) return d;
    return polscamul(a,d,invmod(a[d],p),p);
}

void polcopy( LONG *A, LONG d, LONG *B )
{   int i;
    for( i=0; i<=d; i++) B[i]=A[i];
    return;
}

LONG  polgcd( LONG *A, LONG n, LONG *B, LONG m, LONG p )
// compute the G = gcd(A(x),B(x)) mod p and put the output in A
{   LONG dr,*R1,*R2,*W;
    if( n==-1 && m==-1 ) { printf("A and B must be non-zero\n"); exit(1); }
    if( m==-1 ) { polmonic(A,n,p); return n; } // gcd(A,0) = monic(A)
    if( n==-1 ) { polmonic(B,m,p); return m; } // gcd(0,B) = monic(B)
    R1=A; R2=B;
    while (1)
       {  dr= poldiv(R1,n,R2,m,0,p); // printf("dr:=%lld \n",dr);
           if (dr<0) {
              polmonic(R2,m,p);
              if( A!=R2 ) polcopy(R2, m, A);
              return m;
           }
           W = R1; R1 = R2; R2 = W; n=m; m=dr;
       }
}


LONG Interp(LONG *A,LONG *Y,LONG n,LONG *V, LONG p)
{   LONG d,i,k,j,s;
     if (n==0) return -1;
     V[0]=Y[0];
     if (n==1) { if (V[0]==0) return -1; else return 0;}
     for (j=1; j<=n-1; j++)
     {   k=submod(A[j],A[0],p);
         for (s=V[0],i=1; i<=j-1; i++)
         {   s=addmod(mulmod(k, V[i], p), s, p);
             k=mulmod(submod(A[j], A[i],p), k, p);
         }
         if( k==0 ) { printf("x coordinates must be distinct\n"); exit(1); }
         V[j]=mulmod(submod(Y[j],s,p),invmod(k, p),p);
     }
    for (d=n-1;d>=0 && V[d]==0; d--); //Computing degree d.
    for (i=1; i<=d; i++){
        for (j=d-i; j<=d-1; j++) V[j]= submod(V[j], mulmod(A[d-i],V[j+1],p),p);
    }
    return d;
}




 LONG polGCDEX( LONG *A, int n, LONG *B, int m, LONG *N, LONG *DegN, LONG *D, LONG *DegD, LONG *Temp,LONG p )
 {   /* A is the Treemul pol of deg n>0
        B is the Interp poly of deg m<n.
      */
     LONG dr,*R1,*R2,*W, *Q,*T1,*T2,*T3,dT3,dT1,dT2,Qmax,dN,dD,inv,dn;
     if(n<1 ) return 0; //Qmax=0 => algorithm fails.
     R1=A; R2=B;
     T1=Temp; T2=T1+n+1; T3=T2+n+1;// |T1|=|T2|=|T3|<= n+1. |Temp|=4n+4-m Qmax=m? Qmax= max(n-m,m);
     Q=T3+n+1;
     T1[0]=0; T2[0]=1;
     dT1=-1;  dT2=0;
     Qmax=1;
     while (1)
     {      dr= poldiv(R1,n,R2,m,Q,p); // R1=R1 rem R2
            dT3=polmul(T2,dT2,Q,n-m,T3,p); // T3=T2*Q
            dT3=polsub(T1, dT1, T3, dT3, T3, p);// T3=T1-T3=T1-T2*Q;
            if (n-m > Qmax) {
                Qmax=n-m;
                dN=m; dD=dT2; DegN[0]=dN; DegD[0]=dD;
                polcopy(T2, dD, D);
                polcopy(R2, dN, N);
                dn=polgcd(N, dN, D, dD, p); // In our Gcd routine, dn=0 => (N,D)=1.
                if(dn==0){  //printf("(N,D)=1\n");
                            polcopy(T2, dD, D);polcopy(R2, dN, N);
                            inv=invmod(D[0], p);
                            polscamul(N, dN, inv, p); polscamul(D, dD, inv, p); // The poly N and D normalized by LC(Denom)
                }
                else if( dn!=0 || Qmax<=1) printf("(N,D)!=1\n");// => ( No suitable N and D were found)
            }
            if (dr<0) {// if (m==0) printf("(A,B)=1\n");
                       // if(Qmax<=1) printf("FAIL\n");
                        return Qmax;
                      }
            W = R1; R1 = R2; R2 = W; n=m; m=dr;
            W=T1; T1=T2; T2=T3; T3=W; dT1=dT2; dT2=dT3;
            
    }
 }