void conv(vec_ZZ& x, const vec_zz_p& a)
{
   long n = a.length();
   x.SetLength(n);
   long i;
   for (i = 0; i < n; i++)
      x[i] = rep(a[i]);
}

void ProjectPowers(vec_zz_p& x, const vec_zz_p& a, long k,
                   const zz_pX& h, const zz_pXModulus& F)

{
   if (a.length() > F.n || k < 0) Error("ProjectPowers: bad args");

   if (k == 0) {
      x.SetLength(0);
      return;
   }

   long m = SqrRoot(k);

   zz_pXArgument H;

   build(H, h, F, m);
   ProjectPowers(x, a, k, H, F);
}

void InnerProduct(zz_p& x, const vec_zz_p& a, const vec_zz_p& b,
                  long offset)
{
   if (offset < 0) Error("InnerProduct: negative offset");
   if (NTL_OVERFLOW(offset, 1, 0)) Error("InnerProduct: offset too big");

   long n = min(a.length(), b.length()+offset);
   long i;
   zz_p accum, t;

   clear(accum);
   for (i = offset; i < n; i++) {
      mul(t, a[i], b[i-offset]);
      add(accum, accum, t);
   }

   x = accum;
}

void MinPolySeq(zz_pX& h, const vec_zz_p& a, long m)
{
   if (m < 0 || NTL_OVERFLOW(m, 1, 0)) Error("MinPoly: bad args");
   if (a.length() < 2*m) Error("MinPoly: sequence too short");

   if (m > NTL_zz_pX_BERMASS_CROSSOVER)
      GCDMinPolySeq(h, a, m);
   else
      BerlekampMassey(h, a, m);
}

void FastTraceVec(vec_zz_p& S, const zz_pX& f)
{
   long n = deg(f);

   if (n <= 0) 
      Error("FastTraceVec: bad args");

   if (n == 0) {
      S.SetLength(0);
      return;
   }

   if (n == 1) {
      S.SetLength(1);
      set(S[0]);
      return;
   }
   
   long i;
   zz_pX f1;

   f1.rep.SetLength(n-1);
   for (i = 0; i <= n-2; i++)
      f1.rep[i] = f.rep[n-i];
   f1.normalize();

   zz_pX f2;
   f2.rep.SetLength(n-1);
   for (i = 0; i <= n-2; i++)
      mul(f2.rep[i], f.rep[n-1-i], i+1);
   f2.normalize();

   zz_pX f3;
   InvTrunc(f3, f1, n-1);
   MulTrunc(f3, f3, f2, n-1);

   S.SetLength(n);

   S[0] = n;
   for (i = 1; i < n; i++)
      negate(S[i], coeff(f3, i-1));
}

long IsZero(const vec_zz_p& a)
{
   long n = a.length();
   long i;

   for (i = 0; i < n; i++)
      if (!IsZero(a[i]))
         return 0;

   return 1;
}

void InnerProduct(zz_p& x, const vec_zz_p& a, const vec_zz_p& b)
{
   long n = min(a.length(), b.length());
   long i;

   long accum, t;
   long p = zz_p::modulus();
   mulmod_t pinv = zz_p::ModulusInverse();

   const zz_p *ap = a.elts();
   const zz_p *bp = b.elts();

   accum = 0;
   for (i = 0; i < n; i++) {
      t = MulMod(rep(ap[i]), rep(bp[i]), p, pinv);
      accum = AddMod(accum, t, p);
   }

   x.LoopHole() = accum;
}

void cvt(vec_GF2& x, const vec_zz_p& a)
{
    long n = a.length();

    x.SetLength(n);

    long i;

    for (i = 0; i < n; i++)
        x.put(i, rep(a[i]));
}

void VectorCopy(vec_zz_p& x, const vec_zz_p& a, long n)
{
   if (n < 0) LogicError("VectorCopy: negative length");
   if (NTL_OVERFLOW(n, 1, 0)) ResourceError("overflow in VectorCopy");

   long m = min(n, a.length());

   x.SetLength(n);


   const zz_p *ap = a.elts();
   zz_p *xp = x.elts();

  
   long i;

   for (i = 0; i < m; i++)
      xp[i] = ap[i];

   for (i = m; i < n; i++)
      clear(xp[i]);
}

static
void MulAdd(vec_ZZ& x, const ZZ& prod, const vec_zz_p& h)
{
    long n = x.length();
    if (h.length() != n) Error("MulAdd: dimension mismatch");

    ZZ t;
    long i;

    for (i = 0; i < n; i++) {
        mul(t, prod, rep(h[i]));
        add(x[i], x[i], t);
    }
}

void PlainTraceVec(vec_zz_p& S, const zz_pX& ff)
{
   if (deg(ff) <= 0)
      Error("TraceVec: bad args");

   zz_pX f;
   f = ff;

   MakeMonic(f);

   long n = deg(f);

   S.SetLength(n);

   if (n == 0)
      return;

   long k, i;
   zz_p acc, t;

   const zz_p *fp = f.rep.elts();;
   zz_p *sp = S.elts();

   sp[0] = n;

   for (k = 1; k < n; k++) {
      mul(acc, fp[n-k], k);

      for (i = 1; i < k; i++) {
         mul(t, fp[n-i], rep(sp[k-i]));
         add(acc, acc, t);
      }

      negate(sp[k], acc);
   }
}

void add(vec_zz_p& x, const vec_zz_p& a, const vec_zz_p& b)
{
   long n = a.length();
   if (b.length() != n) LogicError("vector add: dimension mismatch");

   long p = zz_p::modulus();

   x.SetLength(n);

   const zz_p *ap = a.elts();
   const zz_p *bp = b.elts();
   zz_p *xp = x.elts();

   long i;
   for (i = 0; i < n; i++)
      xp[i].LoopHole() = AddMod(rep(ap[i]), rep(bp[i]), p);
}

static
void MixedMul(vec_ZZ& x, const vec_zz_p& a, const mat_ZZ& B)
{
    long n = B.NumRows();
    long l = B.NumCols();

    if (n != a.length())
        Error("matrix mul: dimension mismatch");

    x.SetLength(l);

    long i, k;
    ZZ acc, tmp;

    for (i = 1; i <= l; i++) {
        clear(acc);
        for (k = 1; k <= n; k++) {
            mul(tmp, B(k, i), rep(a(k)));
            add(acc, acc, tmp);
        }
        x(i) = acc;
    }
}

void interpolate(zz_pX& f, const vec_zz_p& a, const vec_zz_p& b)
{
   long m = a.length();
   if (b.length() != m) Error("interpolate: vector length mismatch");

   if (m == 0) {
      clear(f);
      return;
   }

   vec_zz_p prod;
   prod = a;

   zz_p t1, t2;

   long k, i;

   vec_zz_p res;
   res.SetLength(m);

   for (k = 0; k < m; k++) {

      const zz_p& aa = a[k];

      set(t1);
      for (i = k-1; i >= 0; i--) {
         mul(t1, t1, aa);
         add(t1, t1, prod[i]);
      }

      clear(t2);
      for (i = k-1; i >= 0; i--) {
         mul(t2, t2, aa);
         add(t2, t2, res[i]);
      }


      inv(t1, t1);
      sub(t2, b[k], t2);
      mul(t1, t1, t2);

      for (i = 0; i < k; i++) {
         mul(t2, prod[i], t1);
         add(res[i], res[i], t2);
      }

      res[k] = t1;

      if (k < m-1) {
         if (k == 0)
            negate(prod[0], prod[0]);
         else {
            negate(t1, a[k]);
            add(prod[k], t1, prod[k-1]);
            for (i = k-1; i >= 1; i--) {
               mul(t2, prod[i], t1);
               add(prod[i], t2, prod[i-1]);
            }
            mul(prod[0], prod[0], t1);
         }
      }
   }

   while (m > 0 && IsZero(res[m-1])) m--;
   res.SetLength(m);
   f.rep = res;
}

void BuildFromRoots(zz_pX& x, const vec_zz_p& a)
{
   long n = a.length();

   if (n == 0) {
      set(x);
      return;
   }

   long k0 = NextPowerOfTwo(NTL_zz_pX_MUL_CROSSOVER)-1;
   long crossover = 1L << k0;

   if (n <= NTL_zz_pX_MUL_CROSSOVER) {
      x.rep.SetMaxLength(n+1);
      x.rep = a;
      IterBuild(&x.rep[0], n);
      x.rep.SetLength(n+1);
      SetCoeff(x, n);
      return;
   }

   long k = NextPowerOfTwo(n);

   long m = 1L << k;
   long i, j;
   long l, width;

   zz_pX b(INIT_SIZE, m+1);

   b.rep = a;
   b.rep.SetLength(m+1);
   for (i = n; i < m; i++)
      clear(b.rep[i]);

   set(b.rep[m]);
   
   fftRep R1(INIT_SIZE, k), R2(INIT_SIZE, k);


   zz_p t1, one;
   set(one);

   vec_zz_p G(INIT_SIZE, crossover), H(INIT_SIZE, crossover);
   zz_p *g = G.elts();
   zz_p *h = H.elts();
   zz_p *tmp;
   
   for (i = 0; i < m; i+= crossover) {
      for (j = 0; j < crossover; j++)
         negate(g[j], b.rep[i+j]);

      if (k0 > 0) {
         for (j = 0; j < crossover; j+=2) {
            mul(t1, g[j], g[j+1]);
            add(g[j+1], g[j], g[j+1]);
            g[j] = t1;
         }
      }
   
      for (l = 1; l < k0; l++) {
         width = 1L << l;

         for (j = 0; j < crossover; j += 2*width)
            mul(&h[j], &g[j], &g[j+width], width);
      
         tmp = g; g = h; h = tmp;
      }

      for (j = 0; j < crossover; j++)
         b.rep[i+j] = g[j];
   }

   for (l = k0; l < k; l++) {
      width = 1L << l;
      for (i = 0; i < m; i += 2*width) {
         t1 = b.rep[i+width];
         set(b.rep[i+width]);
         TofftRep(R1, b, l+1, i, i+width);
         b.rep[i+width] = t1;
         t1 = b.rep[i+2*width];
         set(b.rep[i+2*width]);
         TofftRep(R2, b, l+1, i+width, i+2*width);
         b.rep[i+2*width] = t1;
         mul(R1, R1, R2);
         FromfftRep(&b.rep[i], R1, 0, 2*width-1);
         sub(b.rep[i], b.rep[i], one);
      }
   }

   x.rep.SetLength(n+1);
   long delta = m-n;
   for (i = 0; i <= n; i++)
     x.rep[i] = b.rep[i+delta];

   // no need to normalize
}

void solve(zz_p& d, vec_zz_p& X, 
           const mat_zz_p& A, const vec_zz_p& b)

{
   long n = A.NumRows();

   if (A.NumCols() != n)
      LogicError("solve: nonsquare matrix");


   if (b.length() != n)
      LogicError("solve: dimension mismatch");

   if (n == 0) {
      set(d);
      X.SetLength(0);
      return;
   }

   long i, j, k, pos;
   zz_p t1, t2, t3;
   zz_p *x, *y;

   mat_zz_p M;
   M.SetDims(n, n+1);
   for (i = 0; i < n; i++) {
      for (j = 0; j < n; j++) 
         M[i][j] = A[j][i];
      M[i][n] = b[i];
   }

   zz_p det;
   set(det);

   long p = zz_p::modulus();
   mulmod_t pinv = zz_p::ModulusInverse();

   for (k = 0; k < n; k++) {
      pos = -1;
      for (i = k; i < n; i++) {
         if (!IsZero(M[i][k])) {
            pos = i;
            break;
         }
      }

      if (pos != -1) {
         if (k != pos) {
            swap(M[pos], M[k]);
            negate(det, det);
         }

         mul(det, det, M[k][k]);

         inv(t3, M[k][k]);
         M[k][k] = t3;


         for (i = k+1; i < n; i++) {
            // M[i] = M[i] - M[k]*M[i,k]*t3

            mul(t1, M[i][k], t3);
            negate(t1, t1);

            x = M[i].elts() + (k+1);
            y = M[k].elts() + (k+1);

            long T1 = rep(t1);
            mulmod_precon_t  t1pinv = PrepMulModPrecon(T1, p, pinv); // T1*pinv;
            long T2;

            for (j = k+1; j <= n; j++, x++, y++) {
               // *x = *x + (*y)*t1

               T2 = MulModPrecon(rep(*y), T1, p, t1pinv);
               x->LoopHole() = AddMod(rep(*x), T2, p);
            }
         }
      }
      else {
         clear(d);
         return;
      }
   }

   X.SetLength(n);
   for (i = n-1; i >= 0; i--) {
      clear(t1);
      for (j = i+1; j < n; j++) {
         mul(t2, X[j], M[i][j]);
         add(t1, t1, t2);
      }
      sub(t1, M[i][n], t1);
      mul(X[i], t1, M[i][i]);
   }

   d = det;
}

void mul(vec_zz_p& x, const vec_zz_p& a, const mat_zz_p& B)
{
   long l = a.length();
   long m = B.NumCols();
  
   if (l != B.NumRows())  
      LogicError("matrix mul: dimension mismatch");  

   if (m == 0) { 

      x.SetLength(0);
      
   }
   else if (m == 1) {

      long p = zz_p::modulus();
      mulmod_t pinv = zz_p::ModulusInverse();

      long acc, tmp;
      long k;

      acc = 0;  
      for(k = 1; k <= l; k++) {  
         tmp = MulMod(rep(a(k)), rep(B(k,1)), p, pinv);  
         acc = AddMod(acc, tmp, p);  
      } 

      x.SetLength(1);
      x(1).LoopHole()  = acc;
          
   }
   else {  // m > 1.  precondition


      long p = zz_p::modulus();
      mulmod_t pinv = zz_p::ModulusInverse();

      vec_long::Watcher watch_mul_aux_vec(mul_aux_vec);
      mul_aux_vec.SetLength(m);
      long *acc = mul_aux_vec.elts();

      long j, k;


      const zz_p* ap = a.elts();

      for (j = 0; j < m; j++) acc[j] = 0;

      for (k = 0;  k < l; k++) {
         long aa = rep(ap[k]);
         if (aa != 0) {
            const zz_p* bp = B[k].elts();
            long T1;
            mulmod_precon_t aapinv = PrepMulModPrecon(aa, p, pinv);

            for (j = 0; j < m; j++) {
               T1 = MulModPrecon(rep(bp[j]), aa, p, aapinv);
               acc[j] = AddMod(acc[j], T1, p);
            }
         } 
      }

      x.SetLength(m);
      zz_p *xp = x.elts();
      for (j = 0; j < m; j++)
         xp[j].LoopHole() = acc[j];    
   }
}