// This file is Cilk versions of the code in fftutil8.c
// I copied my fftutil8.c code and parallelized selected routines.
// Copyright, Michael Monagan, March 2020.

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

// Switch from int i,j,n; to LONG i,j,n; to allow FFTs of size 2^31 and 2^32


#define LONG long long int
#define DEBUG 0

/********************************************************************************/
/*  Zp utilities                                                                */
/********************************************************************************/

#define ADD64s fftadd64scilk
#define ADD32s fftadd32scilk
#define SUB64s fftsub64scilk
#define NEG64s fftneg64scilk
#define MUL64s fftmul64scilk
#define MOD64s fftmod64scilk
#define MIN64s fftmin64scilk
#define INV64s modinv64s

LONG ADD64s(LONG a, LONG b, LONG p) { LONG t; t = (a-p)+b; t += (t>>63) & p; return t; }
LONG ADD32s(int a, int b, int p) { int t; t = (a-p)+b; t += (t>>31) & p; return t; }
LONG SUB64s(LONG a, LONG b, LONG p) { LONG t; t = a-b; t += (t>>63) & p; return t; }
LONG NEG64s(LONG a, LONG p) { LONG t; t = -a; t += (t>>63) & p; return t; }
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;
}

#define UINT64 unsigned long long

typedef struct {
  UINT64 s;	/* shift */
  UINT64 v;	/* reciprocal of d */
  UINT64 d0;	/* divisor shifted up */
  UINT64 d1;
} recint;
recint recip1(UINT64  p);
UINT64 mulrec64(UINT64  a, UINT64  b, recint  v);


LONG MOD64s(LONG a,LONG p) { return(a % p); }
LONG MIN64s(LONG a, LONG b) { if(a<b) return a; else return b; }

LONG modinv64s( LONG c, LONG p );
LONG powmod64s( LONG a, LONG n, LONG p );


/********************************************************************************/
/*  Array utilities                                                             */
/********************************************************************************/

LONG * array64s( LONG n );
void vecprint64s( LONG *A, int n );
void VECCOPY64s( LONG *A, LONG *B, LONG n ); // B[i] = A[i]
void VECFILL64s( LONG *A, LONG n, LONG x ); // A[i] = x
void VECSCAMUL64s( LONG *A, LONG n, LONG x, recint P ); // A[i] = x A[i] mod p
void VECZIPMUL64s( LONG *A, LONG *B, LONG n, recint P ); // A[i] = A[i] B[i] mod p
void VECZIP2MUL64s( LONG *A, LONG *B, LONG n, LONG *C, recint P ); // C[i] = A[i] B[i] mod p
void VECPOWMUL64s( LONG *A, LONG n, LONG omega, recint P ); // A[i] = omega^i mod p

#define VZSIZE 262144
cilk LONG VECSCAMUL64scilk( LONG *A, LONG n, LONG x, recint P ) {
   //VECSCAMUL64s(A,n,x,P);
   //return 1;
   LONG i,m,r,dummy;
   if( n<=VZSIZE ) { VECSCAMUL64s(A,n,x,P); return 1; }
   m = n/VZSIZE; r = n-m*VZSIZE;
   for( i=0; i<m; i++ ) { dummy = spawn VECSCAMUL64scilk(A,VZSIZE,x,P); A += VZSIZE; }
   VECSCAMUL64s(A,r,x,P);
   sync;
   return 1;
}
cilk LONG VECZIPMUL64scilk( LONG *A, LONG *B, LONG n, recint P ) {
   LONG i,m,r,dummy;
   //VECZIPMUL64s(A,B,n,P);
   //return 1;
   if( n<=VZSIZE ) { VECZIPMUL64s(A,B,n,P); return 1; }
   m = n/VZSIZE; r = n-m*VZSIZE;
   for( i=0; i<m; i++ ) { dummy = spawn VECZIPMUL64scilk(A,B,VZSIZE,P); A += VZSIZE; B += VZSIZE; }
   VECZIPMUL64s(A,B,r,P);
   sync;
   return 1;
}

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

int poladd64s(LONG *a, LONG *b, LONG *c, int da, int db, LONG p);
int polsqr64s( LONG * A, LONG * C, int d, 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 );
void polprint64s( LONG *A, int d );
void polcopy64s( LONG *A, int d, LONG *B );
LONG poleval64s(LONG *a, int d, LONG x, LONG p);


/********************************************************************************/
/*  FFT utilities and routines                                                  */
/********************************************************************************/

int vecequal64s( LONG *A, LONG *B, int n );

LONG getprimelem64s( LONG p ); // in primitive.c
LONG getomega64s( LONG p, LONG n ); // p is assumed prime > 2

#define FFTCUTOFF 200000
void FFT64s1( LONG n, LONG *a, LONG *W, LONG p, recint P );
void FFT64s2( LONG n, LONG *a, LONG *W, LONG p, recint P );
void FFTwork1( LONG n2, LONG *a, LONG *b, LONG *W, LONG p, recint P );
void FFTwork2( LONG n2, LONG *a, LONG *b, LONG *W, LONG p, recint P );

#define FFTSIZE 262144
cilk LONG FFT64s1workcilk( LONG n2, LONG *a, LONG *b, LONG *W, LONG p, recint P ) {
   LONG i,dummy,m;
   if( n2<=FFTSIZE ) { FFTwork1(n2,a,b,W,p,P); return 1; }
   m = n2/FFTSIZE;
   if( n2-m*FFTSIZE != 0 ) printf("FFTSIZE bug\n");
   for( i=0; i<m; i++ ) {
       dummy = spawn FFT64s1workcilk( FFTSIZE, a, b, W, p, P );
       a += FFTSIZE; b += FFTSIZE; W += FFTSIZE;
   }
   sync;
   return 1;
}

cilk LONG FFT64s1cilk( LONG n, LONG *a, LONG *W, LONG p, recint P ) {
   LONG i,n2,*b;
   if( n<=FFTCUTOFF ) { FFT64s1(n,a,W,p,P); return 1; }
   n2 = n/2;
   b = a+n2;
   //for( i=0; i<n2; i++ ) {
   //   t = SUB64s(a[i],b[i],p);
   //   t = mulrec64(t,W[i],P);
   //   a[i] = ADD64s(a[i],b[i],p);
   //   b[i] = t;
   //}
   //FFTwork1(n2,a,b,W,p,P);
   i = spawn FFT64s1workcilk(n2,a,b,W,p,P);
   sync;
   i = spawn FFT64s1cilk(n2,a,W+n2,p,P);
   i = spawn FFT64s1cilk(n2,b,W+n2,p,P);
   sync;
   return 1;
}


cilk LONG FFT64s2workcilk( LONG n2, LONG *a, LONG *b, LONG *W, LONG p, recint P ) {
   LONG i,dummy,m;
   if( n2<=FFTSIZE ) { FFTwork2(n2,a,b,W,p,P); return 1; }
   m = n2/FFTSIZE;
   if( n2-m*FFTSIZE != 0 ) printf("FFTSIZE bug\n");
   for( i=0; i<m; i++ ) {
       dummy = spawn FFT64s2workcilk( FFTSIZE, a, b, W, p, P );
       a += FFTSIZE; b += FFTSIZE; W += FFTSIZE;
   }
   sync;
   return 1;
}


cilk LONG FFT64s2cilk( LONG n, LONG *a, LONG *W, LONG p, recint P )
{ LONG i,n2,*b;
  if( n<FFTCUTOFF ) { FFT64s2(n,a,W,p,P); return 1; }
  n2 = n/2;
  b = a+n2;
  i = spawn FFT64s2cilk(n2,a,W+n2,p,P);
  i = spawn FFT64s2cilk(n2,b,W+n2,p,P);
  sync;
  //for( i=0; i<n2; i++ ) {
  //    s = a[i];
  //    t = mulrec64(b[i],W[i],P);
  //    a[i] = ADD64s(s,t,p);
  //    b[i] = SUB64s(s,t,p);
  //}
  //FFTwork2(n2,a,b,W,p,P);
  i = spawn FFT64s2workcilk(n2,a,b,W,p,P);
  sync;
  return 1;
}


void FFT64s2even( LONG n, LONG *a, LONG *W, LONG p, recint P );
void FFT64s2even2( LONG n, LONG *a, LONG *W, LONG p, recint P );
void FFT64s2T( LONG n, LONG *a, LONG *W, LONG *T, LONG p, recint P );
void FFTPermute64s( LONG n, LONG *A, LONG *T );
void polmod64s( LONG d, LONG n, LONG *A, LONG *B, LONG p );
void polpad64s( LONG d, LONG n, LONG *A, LONG *B );
void MakeW64( LONG n, LONG w, LONG *W, LONG p );
void MakeWinv64( LONG n, LONG *W, LONG p );


/********************************************************************************/
/*  Fast polynomial multiplication and division                                 */
/********************************************************************************/

int FFTmul64s( LONG * a, LONG * b, LONG * c, int da, int db, LONG *T, LONG p );
cilk int FFTmul64scilk( LONG * a, LONG * b, LONG * c, int da, int db, LONG *T, LONG p )
{  // T has to be length at least 3n where n=2^k > da+db
   int d;  LONG i,n,dc;  LONG dummy,w,t,*A,*B,*W;  recint P;

   if( da==-1 || db==-1 ) return -1;
   P = recip1(p);
   dc = da+db;
   for( n=1; n<dc; n*=2 );    // printf("da=%d db=%d n=%d\n",da,db,n);
   if( n<FFTCUTOFF ) return FFTmul64s(a,b,c,da,db,T,p);
   if( n==dc ) { // avoid an FFT with size 2 dc
       d = spawn FFTmul64scilk( a, b+1, c+1, da, db-1, T, p );
       sync;
       c[0] = mulrec64(a[0],b[0],P);
       for( i=1; i<=da; i++ ) c[i] = ADD64s(c[i],mulrec64(b[0],a[i],P),p);
       for( dc=da+db; dc>=0 && c[dc]==0; dc-- );
       return dc;
   }
   w = getomega64s(p,n);                   //printf("w := %lld;\n",w);
   if( w==0 ) { printf("omega does not exist: dc=%d  n=%d\n",dc,n); exit(1); }
   A = T; B = T+n; W = B+n;
   MakeW64(n,w,W,p);                            //printf("W := "); vecprint32s(W,n);
   polpad64s(da,n,A,a);                         //printf("A := "); vecprint32s(A,n);
   dummy = spawn FFT64s1cilk(n,A,W,p,P);        //printf("F := "); vecprint32s(A,n);
   if( a==b ) B = A; // c = a^2
   else {
       polpad64s(db,n,B,b);                     //printf("B := "); vecprint32s(B,n);
       dummy = spawn FFT64s1cilk(n,B,W,p,P);    //printf("F := "); vecprint32s(B,n);
   }
   sync;
   for( i=0; i<n; i++ ) A[i] = mulrec64(A[i],B[i],P);  //printf("AB := "); vecprint32s(A,n);
   MakeWinv64(n,W,p);                           //printf("W := "); vecprint32s(W,n);
   dummy = spawn FFT64s2cilk(n,A,W,p,P);        //printf("C := "); vecprint32s(A,n);
   sync;
   t = INV64s(n,p);
   for( i=0; i<=dc; i++ ) c[i] = mulrec64(t,A[i],P);
   while( dc>=0 && c[dc]==0 ) dc--;
   return dc;
}


int FFTpolmul64s( LONG *a, LONG *b, LONG *c, LONG da, LONG db, LONG p );
cilk int FFTpolmul64scilk( LONG *a, LONG *b, LONG *c, LONG da, LONG db, LONG p )
{
    int dc; LONG n,*T;
    if( da<0 || db<0 ) return -1;
    if( da<20 || db<20 || (LONG) da*db < 4096 ) return polmul64s(a,b,c,da,db,p);
    dc = da+db;
    for( n=1; n<=dc; n*=2 ); // FFT must hold a x b which has degree dc so n must be > dc.
    if( n<FFTCUTOFF ) return FFTpolmul64s(a,b,c,da,db,p);
    T = array64s(3*n);
    dc = spawn FFTmul64scilk( a, b, c, da, db, T, p );
    sync;
    free(T);
    return dc;
}


/*********************   Joris' root finding algorithm    *********************/


void FFTgraeffe64s( LONG * f, LONG * g, LONG d, LONG *T, LONG p );
void tangGraeffe64s( LONG *f, LONG *g, LONG d, LONG p );
void compress64s( LONG *a, LONG n, LONG p );
void FFTtangGraeffeWork( LONG n2, LONG *A, LONG *C, LONG p, recint P );

void FFTtangGraeffe64s( LONG *f, LONG *g, LONG d, int k, LONG *T, LONG p );
cilk int FFTtangGraeffe64scilk( LONG *f, LONG *g, LONG d, int k, LONG *T, LONG p )
{
   // This version composes the Tangent Graeffe transform k times
   //
   // Given h(x) = f(x) + e g(x) where deg(f)=d and deg(g)=d-1
   // Compute r(x,e) where r(x^2,e) = h(x) h(-x) mod e^2
   // Via P(x^2) = f(x) f(-x) and Q(x^2) = g(x) f(-x) + f(x) g(-x)
   // I decided to overwrite f with P(x) and g with Q(x) so that the output is f(x) + g(x) e
   // The working storage vector T needs to be of size 3n where n=2^k>2d
   // There is an optimization available for the inverse transform of vectors
   // e.g. (n=8) [a0, a0, a1, a1, a2, a2, a3, a3] 
   LONG j;
   LONG n,n2,n4;
   LONG dummy, w, t, *A, *C, *W; 
   recint P;
   P = recip1(p);
   if( d<1 ) { printf("degree d must be at least 1\n"); exit(1); }
   for( n=1; n<=2*d; n*=2 );                        //printf("n := %d;\n",n);
   if( n<FFTCUTOFF ) { FFTtangGraeffe64s(f,g,d,k,T,p); return 1; }
   w = getomega64s(p,n);                            //printf("w := %lld;\n",w);
   if( w==0 ) { printf("omega does not exist: d=%d  n=%d\n",d,n); exit(1); }
   n2 = n/2;
   n4 = n/4;
   A = T; C = T+n; W = C+n; // B = W+n; D = B+n;
   MakeW64(n,w,W,p);                                //printf("W := "); vecprint64s(W,n);
   polpad64s(d,n,A,f);                              //printf("A := "); vecprint64s(A,n);
   polpad64s(d-1,n,C,g);                            //printf("C := "); vecprint64s(C,n-1);
   dummy = spawn FFT64s1cilk(n,A,W,p,P); // A = FFT1(f(x))  //printf("A := "); vecprint64s(A,n);
   dummy = spawn FFT64s1cilk(n,C,W,p,P); // C = FFT1(g(x))  //printf("C := "); vecprint64s(A,n);
   sync;
for( j=0; j<k; j++ ) {
   //if( n>1000000 ) printf(" TG: j=%d\n", j );
   if( j>0 ) { // double the order
       MakeWinv64(n,W,p); // reset it
       //VECZIPMUL64s(A+n2,W,n2,P); //for( i=0; i<n2; i++ ) A[n2+i] = mulrec64(A[n2+i],W[i],P);
       //VECZIPMUL64s(C+n2,W,n2,P); //for( i=0; i<n2; i++ ) C[n2+i] = mulrec64(C[n2+i],W[i],P);
       dummy = spawn VECZIPMUL64scilk(A+n2,W,n2,P);
       dummy = spawn VECZIPMUL64scilk(C+n2,W,n2,P);
       sync;
       dummy = spawn FFT64s1cilk(n2,A+n2,W+n2,p,P);
       dummy = spawn FFT64s1cilk(n2,C+n2,W+n2,p,P);
       sync;
   }
   // The FFT permutation for [0 1 2 3 4 5 6 7] is [0 4 2 6 1 5 3 7]
   // The FFT permutation for f(-x) is [4 0 6 2 5 1 7 3] so just interchange
   // for( i=0; i<n; i+=2 ) { B[i] = A[i+1]; B[i+1] = A[i]; } // B = FFT1(f(-x))
   // for( i=0; i<n; i+=2 ) { D[i] = C[i+1]; D[i+1] = C[i]; } // D = FFT1(g(-x))
   // But we don't need to explicitly do that because we know where the negatives are
   //for( i=0; i<n2; i++ ) // C = BC + AD = FFT( f(-x) g(x) + f(x) g(-x) )
   //     //C[i] = ADD64s( mulrec64(C[2*i],B[2*i],P), mulrec64(D[2*i],A[2*i],P), p );
   //     C[i] = ADD64s( mulrec64(C[2*i],A[2*i+1],P), mulrec64(C[2*i+1],A[2*i],P), p );
   //for( i=0; i<n2; i++ ) // A = A B = FFT( f(x) f(-x) )
   //     //A[i] = mulrec64(A[2*i],B[2*i],P); 
   //     A[i] = mulrec64(A[2*i],A[2*i+1],P); 
   FFTtangGraeffeWork( n2, A, C, p, P );
   // A = [a0,--,a1,--,a2,--,a3,--] and C = [c0,--,c1,--,c2,--,c3,--]
   VECCOPY64s(A,A+n2,n2);  // A = [a0,a1,a2,a3,a0,a1,a2,a3]
   VECCOPY64s(C,C+n2,n2);  // C = [c0,c1,c2,c3,c0,c1,c2,c3]
   // Inverse FFT on the copies
   MakeWinv64(n,W,p);
   dummy = spawn FFT64s2cilk(n2,A+n2,W+n2,p,P);
   dummy = spawn FFT64s2cilk(n2,C+n2,W+n2,p,P);
   sync;
   t = INV64s(n2,p);
   if( d&1 ) t = NEG64s(t,p); // to make lc(f) = +1 and lc(g) = +1
   //VECSCAMUL64s(A+n2,n2,t,P);
   //VECSCAMUL64s(C+n2,n2,t,P);
   dummy = spawn VECSCAMUL64scilk(A+n2,n2,t,P);
   dummy = spawn VECSCAMUL64scilk(C+n2,n2,t,P);
   sync;
}
   VECCOPY64s(A+n2,f,d+1);
   VECCOPY64s(C+n2,g,d);
   return 1;
}


LONG bluestein64s( LONG *a, LONG d, LONG s, LONG *v, LONG p );
LONG FFTbluestein64s( LONG *a, LONG d, LONG s, LONG *v, LONG *T, LONG p );
cilk LONG FFTbluestein64scilk( LONG *a, LONG d, LONG s, LONG *v, LONG *T, LONG p )
{
   // a(x) is a polynomial of degree d
   // compute v = [ a(w^i) mod p, 0 <= i < s ] for w^s=1 via Bluestein
   // BUG:  the computation (i*i)%(2s) can overflow in an int type
   // MBM:  the (i*i)%(2s) can be replaced using i*i = (i-1)*(i-1)+2i-1

   LONG i,j,k,N;
   LONG w,omega,t,*A,*B,*W; 
   recint P;

   if( s<2*FFTCUTOFF ) return FFTbluestein64s(a,d,s,v,T,p);
   P = recip1(p);
   //printf("Bluestein: d=%d  s=%d\n",d,s);
   if( (p-1)%s != 0 ) { printf("s must divide p-1\n"); exit(1); }
   w = getomega64s(p,2*s);                          //printf("w := %lld;\n",w);
   if( w==0 ) { printf("omega does not exist: s=%d\n",s); exit(1); }
   omega = mulrec64(w,w,P);
   for( N=2; N<2*s; N=2*N );                        // printf("Bluestein: N = %d\n",N);
   if( (p-1)%N != 0 ) { printf("N must divide p-1\n"); exit(1); }
   // T = array64s(3*N);  T must be at least this big
   W = T; A = T+N; B = A+N;
   //W[0] = 1; W[s] = p-1;
   //for( i=1; i<s; i++ ) { W[i] = mulrec64(w,W[i-1],P); W[s+i] = NEG64s(W[i],p); }
   VECPOWMUL64s(W,s,w,P); for( i=0; i<s; i++ ) W[s+i] = NEG64s(W[i],p);
   v[0] = 1;
   //for( i=1; i<s; i++ ) v[i] = W[((LONG) i*i)%(2*s)];      //printf("v := "); vecprint64s(v,s);
   for( k=0,i=1,j=1; i<s; i++,j+=2 ) { k = ADD32s(j,k,2*s); v[i] = W[k]; }
   //for( i=0; i<=d; i++ ) A[i] = mulrec64(a[i],v[i],P);
   VECZIP2MUL64s(a,v,d+1,A,P);
   for( i=d+1; i<N; i++ ) A[i] = 0;                 //printf("A := "); vecprint64s(A,N);
   // compute W = [1,w^(-1),w^(-2),...,w^(1-2s)] from W = [1,w,w^2,...,w^(2s-1)]
   for( i=1,j=2*s-1; i<s; i++,j-- ) { t = W[i]; W[i] = W[j]; W[j] = t; }
   B[0] = 1;
   //for( i=1; i<s; i++ ) B[i] = W[((LONG) i*i)%(2*s)];
   for( k=0,i=1,j=1; i<s; i++,j+=2 ) { k = ADD32s(j,k,2*s); B[i] = W[k]; }
   for( i=s; i<N; i++ ) B[i] = 0;
   for( i=1; i<s; i++ ) B[N-i] = B[i];              //printf("B := "); vecprint64s(B,N);
   // We will multiply A B mod x^N - 1 using an FFT of order N
   w = getomega64s(p,N);                            //printf("w := %lld;\n",w);
   if( w==0 ) { printf("omega does not exist: N=%d\n",N); exit(1); }
   MakeW64(N,w,W,p);
   i = spawn FFT64s1cilk(N,A,W,p,P);                //printf("FA := "); vecprint64s(A,N);
   i = spawn FFT64s1cilk(N,B,W,p,P);                //printf("FB := "); vecprint64s(B,N);
   sync;
   VECZIPMUL64s(A,B,N,P); // for( i=0; i<N; i++ ) A[i] = mulrec64(A[i],B[i],P);
   MakeWinv64(N,W,p);
   i = spawn FFT64s2cilk(N,A,W,p,P);
   sync;
   t = INV64s(N,p);
   //for( i=0; i<s; i++ ) { w = mulrec64(t,A[i],P); v[i] = mulrec64(v[i],w,P); }
   VECSCAMUL64s(A,s,t,P); // A[i] = t A[i] mod p
   VECZIPMUL64s(v,A,s,P); // v[i] = v[i] A[i] mod p
   return omega; // omega^s = 1
}


/*****************************    Fast division    ******************************/


void polinv64s( LONG *f, int n, LONG *y, LONG *T, LONG p );

#define INVCUTOFF 128
void FFTpolinvmod64s( LONG *f, LONG d, LONG n, LONG *y, LONG *W, LONG *Winv, LONG *T, LONG p );
cilk int FFTpolinvmod64scilk( LONG *f, LONG d, LONG n, LONG *y, LONG *W, LONG *Winv, LONG *T, LONG p )
{
    LONG i,m,n2;
    recint P;
    LONG *Y,*F,*M,ni;
    if( n<INVCUTOFF ) { polinv64s( f, MIN64s(n,d), y, T, p ); return 1; }
    Y = T;
    F = T+n;
    n2 = n/2;
    P = recip1(p);
    FFTpolinvmod64s( f, d, n2, y, W+n2, Winv+n2, T, p );
    VECCOPY64s(y,Y,n2); VECFILL64s(Y+n2,n2,0);
    FFT64s1(n,Y,W,p,P);
    m = MIN64s(d,n); VECCOPY64s(f,F,m); VECFILL64s(F+m,n-m,0);
    FFT64s1(n,F,W,p,P);
    VECZIPMUL64s(F,Y,n,P); //for( i=0; i<n; i++ ) F[i] = MUL(F[i],Y[i],p);
    FFT64s2(n,F,Winv,p,P);
    ni = modinv64s(n,p);
    VECSCAMUL64s(F,n,ni,P);
    // M = "middle product" is in second half of F
    M = F+n2; VECCOPY64s(M,F,n2); VECFILL64s(M,n2,0);
    FFT64s1(n,F,W,p,P);
    VECZIPMUL64s(F,Y,n,P); // for( i=0; i<n; i++ ) F[i] = MUL(F[i],Y[i],p);
    FFT64s2(n,F,Winv,p,P);
    VECSCAMUL64s(F,n,ni,P);
    for( i=0; i<n2; i++ ) y[n2+i] = NEG64s(F[i],p);      // yk = yk + x^n2 F
    return 1;
}

void FFTinv64s( LONG * f, LONG m, LONG * y, LONG *T, LONG p );
cilk int FFTinv64scilk( LONG * f, LONG m, LONG * y, LONG *T, LONG p )
{  // y must be of length m
   // T must be of length 5n where n=2^k and n >= m
   LONG i,n;
   LONG w,winv,*W,*Winv,*Y;
   if( m==1 ) { y[0] = INV64s(f[0],p); return 1; }
   for( n=1; n<m; n*=2 );                       //printf("n := %d;\n",n);
   w = getomega64s(p,n);                        //printf("w := %lld;\n",w);
   if( w==0 ) { printf("omega does not exist  n=%d\n",n); exit(1); }
   W = T;
   MakeW64(n,w,W,p);
   winv = modinv64s(w,p);
   Winv = T + n; 
   MakeW64(n,winv,Winv,p);
   Y = T + 2*n;
   FFTpolinvmod64s( f, m, n, Y, W, Winv, T+3*n, p );
   for( i=0; i<m; i++ ) y[i] = Y[i];
   return 1;
}

int FFTpoldivinp64s( LONG *a, LONG *b, LONG da, LONG db, LONG *Q, LONG *T, LONG p );
cilk int FFTpoldivinp64scilk( LONG *a, LONG *b, LONG da, LONG db, LONG *Q, LONG *T, LONG p ) {
// Inplace a div b using the FFT
// Q must be an array of size 4(dq+1) and T of size 4n where n=2^k>da 
    LONG i,dr,dq,n;
    LONG *ra, *rb, *rq, *y, *q, *r;
    if( db<0 ) { printf("division by zero\n"); exit(1); }
    if( da<db ) return da;
    dq = da-db;
    if( db<16 || dq<16 || (LONG) db*dq < 40000 ) return poldiv64s(a,b,da,db,p);
    q = a + db; // this is where q will go
    for( n=1; n<=da; n*=2 ); // FFT must hold q x b which has degree da so n must be > da.
    //Q = array64s(4*dq+4);
    //T = array64s(4*n);
    ra = Q; y = Q + dq + 1; rb = Q + 2*dq + 2;
    for( i=0; i<=dq; i++ ) ra[i] = a[da-i]; // ra = recip(a)
//printf("p:="); printf("%lld;\n",p);
//printf("ra:="); polprint64s(ra,dq);
    for( i=0; i<=MIN64s(db,dq); i++ ) rb[i] = b[db-i]; // rb = recip(b)
//printf("rb:="); polprint64s(rb,dq);
    while( i<=dq ) rb[i++] = 0;
    FFTinv64s( rb, dq+1, y, T, p ); // compute 1/rb to O(x^(dq+1))
//printf("inv:="); polprint64s(y,dq);
    rq = Q + 2*dq + 2;
    FFTmul64s( ra, y, rq, dq, dq, T, p ); // T must be size 3n
//printf("rq:="); polprint64s(rq,2*dq);
    for( i=0; i<=dq; i++ ) q[i] = rq[dq-i];
    r = T;
    FFTmul64s( b, q, r, db, dq, T+n, p ); // T must be size 4n
    for( i=0; i<db; i++ ) a[i] = SUB64s(a[i],r[i],p); // copy the remainder into a
    for( dr=db-1; dr>=0 && a[dr]==0; dr-- ); // compute deg(r)
    return dr;
}
    

int FFTpoldiv64s( LONG *a, LONG *b, LONG da, LONG db, LONG p );
cilk int FFTpoldiv64scilk( LONG *a, LONG *b, LONG da, LONG db, LONG p ) {
    LONG dr,dq,n;
    LONG *Q, *T;
    if( db<0 ) { printf("division by zero\n"); exit(1); }
    if( da<db ) return da;  
    dq = da-db; 
    if( db<16 || dq<16 || (LONG) db*dq < 40000 ) return poldiv64s(a,b,da,db,p);
    for( n=1; n<=da; n*=2 ); // FFT must hold q x b which has degree da so n must be > da.
//printf("FFTpoldiv64s: da=%d  n=%d  p=%lld\n",da,n,p);
    Q = array64s(4*dq+4);
    T = array64s(4*n);
    dr = FFTpoldivinp64s( a, b, da, db, Q, T, p );
//printf("FFTpoldiv64s: dr=%d\n",dr);
    free(Q);
    free(T);
    return dr;
}


/*****************************  TREE MUL ALGORITHM  ********************************/


LONG treemulspace( LONG n, LONG m );

#define CUTOFF 16
#define LAMBDACUTOFF 100000
void fastLambdarec( LONG *alpha, LONG n, LONG *f, LONG *T, LONG s, LONG *W, LONG p );
cilk int fastLambdareccilk( LONG *alpha, LONG n, LONG *f, LONG *T, LONG s, LONG *W, LONG Wsize, LONG p ) {
    int dummy;
    LONG m;
    if( n<LAMBDACUTOFF ) { fastLambdarec(alpha,n,f,T,s,W,p); return 1; }
    m = n/2;
    dummy = spawn fastLambdareccilk( alpha, m, T, f, s/2, W, Wsize/2, p );
    dummy = spawn fastLambdareccilk( alpha+m, n-m, T+s/2, f+s/2, s/2, W+Wsize/2, Wsize/2, p );
    sync;
    // FFTpolmul64s( T, T+s/2, f, m, n-m, p );
    if( n<200 ) { polmul64s( T, T+s/2, f, m, n-m, p ); return 1; }
    //FFTmul64s( T, T+s/2, f, m, n-m, W, p );
    dummy = spawn FFTmul64scilk( T, T+s/2, f, m, n-m, W, p );
    sync;
    return 1;
}

void fastLambda( LONG * v, LONG n, LONG *f, LONG p );
cilk int fastLambdacilk( LONG * v, LONG n, LONG *f, LONG p ) {
// v = [a1,a2,a3,...,an] compute lambda(x) = (x-a1)(x-a2)...(x-an)
    int dummy;
    LONG s, N, Wsize;
    LONG *T, *F, *W;
    if( n<LAMBDACUTOFF ) { fastLambda(v,n,f,p); return 1; }
    s = treemulspace( n, CUTOFF ); //printf("s = %d\n",s);
    T = array64s(s);
    F = array64s(s);
    for( N=1; N<=n; N*=2 ); // FFT must hold lambda(x) of degree n so N > n.
    //printf("N = %d\n",N);
    Wsize = 3*N;
    W = array64s(Wsize);
    dummy = spawn fastLambdareccilk( v, n, F, T, s, W, Wsize, p );
    sync;
    polcopy64s(F,n,f);
    free(F);
    free(W);
    free(T);
    return 1;
}


/*****************************  Fast change of basis  ******************************/


void makedivisors( LONG alpha, LONG n, LONG *T, LONG p );
void makepowers( LONG alpha, LONG n, LONG *T, LONG p );
void recbase( LONG *f, LONG n, LONG alpha, LONG *T, LONG m, LONG p );
void recbase2( LONG *f, LONG n, LONG alpha, LONG *T, LONG m, LONG *W, LONG *S, LONG p );
int divisorsspace( LONG n );
void changebase( LONG *f, LONG n, LONG alpha, LONG p );
void changebase2( LONG *f, LONG n, LONG alpha, LONG p );
void polrev64s( LONG *f, LONG d );


void changebase3( LONG *f, LONG n, LONG alpha, LONG p );
cilk LONG changebase3cilk( LONG *f, LONG n, LONG alpha, LONG p ) {
//  Input f a polynomial of degree n in Zp[x]
//  Output f(x+alpha) mod p in the array f
//  This version computes f(x+alpha) in O( M(n) ) instead of O( M(n) log n )
//  Since the method computes 1/i! mod p we cannot use it for p<=n.
//  Computing inverses is expensive so some thought is needed for that.
//  Joris suggested precomputing (1/i!)^(-1) and going backwards.
//  Since the we have to compute i! first anyway this is free.
    LONG i,fac,inv,*g,*h;
    int dummy;
    recint P;
    if( p<=n || n<3000 ) { changebase2(f,n,alpha,p); return 1; }
    if( n<FFTCUTOFF ) { changebase3(f,n,alpha,p); return 1; }
    P = recip1(p);
    fac = 1;
    for( i=1; i<=n; i++ ) {
        // compute f[i] = f[i]*i! mod p
        fac = mulrec64(i,fac,P);
        f[i] = mulrec64(fac,f[i],P);
    }
    inv = modinv64s(fac,p); // = 1/n! mod p
    polrev64s(f,n);
    g = array64s(n+1);
    g[0] = 1;
    for( i=1; i<=n; i++ ) {
        // compute g[i] = alpha^i mod p
        g[i] = mulrec64(alpha,g[i-1],P);
    }
    fac = inv;
    for( i=n; i>0; i-- ) {
        // compute g[i] = alpha^i/i! mod p
        g[i] = mulrec64(fac,g[i],P);
        fac = mulrec64((LONG) i,fac,P); 
    }
    h = array64s(2*n+1);
    dummy = spawn FFTpolmul64scilk(f,g,h,n,n,p);
    sync;
    polrev64s(h,n);
    f[0] = h[0];
    fac = inv;
    for( i=n; i>0; i-- ) {
        // compute f[i] = h[i]/i! mod p
        f[i] = mulrec64(fac,h[i],P);
        fac = mulrec64(fac,(LONG) i,P);
    }
    free(g);
    free(h);
    return 1;
}