add(a[i], t, a[i-1]);      }      mul(a[0], a[0], b);   }} void mul(ZZ_p* x, const ZZ_p* a, const ZZ_p* b, long n){   static ZZ t, accum;   long i, j, jmin, jmax;   long d = 2*n-1;   for (i = 0; i <= d; i++) {      jmin = max(0, i-(n-1));      jmax = min(n-1, i);      clear(accum);      for (j = jmin; j <= jmax; j++) {         mul(t, rep(a[j]), rep(b[i-j]));         add(accum, accum, t);      }      if (i >= n) {         add(accum, accum, rep(a[i-n]));         add(accum, accum, rep(b[i-n]));      }      conv(x[i], accum);   }}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_FFT_CROSSOVER);   long crossover = 1L << k0;   if (n <= 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 eval(ZZ_p& b, const ZZ_pX& f, const ZZ_p& a)// does a Horner evaluation{   ZZ_p acc;   long i;   clear(acc);   for (i = deg(f); i >= 0; i--) {      mul(acc, acc, a);      add(acc, acc, f.rep[i]);   }   b = acc;}void eval(vec_ZZ_p& b, const ZZ_pX& f, const vec_ZZ_p& a)// naive algorithm:  repeats Horner{   if (&b == &f.rep) {      vec_ZZ_p bb;      eval(bb, f, a);      b = bb;      return;   }   long m = a.length();   b.SetLength(m);   long i;   for (i = 0; i < m; i++)       eval(b[i], f, a[i]);}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;}NTL_vector_impl(ZZ_pX,vec_ZZ_pX)NTL_eq_vector_impl(ZZ_pX,vec_ZZ_pX)NTL_io_vector_impl(ZZ_pX,vec_ZZ_pX)   void InnerProduct(ZZ_pX& x, const vec_ZZ_p& v, long low, long high,                    const vec_ZZ_pX& H, long n, ZZVec& t){   static ZZ s;   long i, j;   for (j = 0; j < n; j++)      clear(t[j]);   high = min(high, v.length()-1);   for (i = low; i <= high; i++) {      const vec_ZZ_p& h = H[i-low].rep;      long m = h.length();      const ZZ& w = rep(v[i]);      for (j = 0; j < m; j++) {         mul(s, w, rep(h[j]));         add(t[j], t[j], s);      }   }   x.rep.SetLength(n);   for (j = 0; j < n; j++)      conv(x.rep[j], t[j]);   x.normalize();}void CompMod(ZZ_pX& x, const ZZ_pX& g, const ZZ_pXArgument& A,              const ZZ_pXModulus& F){   if (deg(g) <= 0) {      x = g;      return;   }   ZZ_pX s, t;   ZZVec scratch(F.n, ZZ_pInfo->ExtendedModulusSize);   long m = A.H.length() - 1;   long l = ((g.rep.length()+m-1)/m) - 1;   ZZ_pXMultiplier M;   build(M, A.H[m], F);   InnerProduct(t, g.rep, l*m, l*m + m - 1, A.H, F.n, scratch);   for (long i = l-1; i >= 0; i--) {      InnerProduct(s, g.rep, i*m, i*m + m - 1, A.H, F.n, scratch);      MulMod(t, t, M, F);      add(t, t, s);   }   x = t;}void build(ZZ_pXArgument& A, const ZZ_pX& h, const ZZ_pXModulus& F, long m){   if (m <= 0 || deg(h) >= F.n) Error("build: bad args");   if (m > F.n) m = F.n;   long i;   if (ZZ_pXArgBound > 0) {      double sz = ZZ_p::storage();      sz = sz*F.n;      sz = sz + NTL_VECTOR_HEADER_SIZE + sizeof(vec_ZZ_p);      sz = sz/1024;      m = min(m, long(ZZ_pXArgBound/sz));      m = max(m, 1);   }   ZZ_pXMultiplier M;   build(M, h, F);   A.H.SetLength(m+1);   set(A.H[0]);   A.H[1] = h;   for (i = 2; i <= m; i++)       MulMod(A.H[i], A.H[i-1], M, F);}long ZZ_pXArgBound = 0;void CompMod(ZZ_pX& x, const ZZ_pX& g, const ZZ_pX& h, const ZZ_pXModulus& F)   // x = g(h) mod f{   long m = SqrRoot(g.rep.length());   if (m == 0) {      clear(x);      return;   }   ZZ_pXArgument A;   build(A, h, F, m);   CompMod(x, g, A, F);}void Comp2Mod(ZZ_pX& x1, ZZ_pX& x2, const ZZ_pX& g1, const ZZ_pX& g2,              const ZZ_pX& h, const ZZ_pXModulus& F){   long m = SqrRoot(g1.rep.length() + g2.rep.length());   if (m == 0) {      clear(x1);      clear(x2);      return;   }   ZZ_pXArgument A;   build(A, h, F, m);   ZZ_pX xx1, xx2;   CompMod(xx1, g1, A, F);   CompMod(xx2, g2, A, F);   x1 = xx1;   x2 = xx2;}void Comp3Mod(ZZ_pX& x1, ZZ_pX& x2, ZZ_pX& x3,               const ZZ_pX& g1, const ZZ_pX& g2, const ZZ_pX& g3,              const ZZ_pX& h, const ZZ_pXModulus& F){   long m = SqrRoot(g1.rep.length() + g2.rep.length() + g3.rep.length());   if (m == 0) {      clear(x1);      clear(x2);      clear(x3);      return;   }   ZZ_pXArgument A;   build(A, h, F, m);   ZZ_pX xx1, xx2, xx3;   CompMod(xx1, g1, A, F);   CompMod(xx2, g2, A, F);   CompMod(xx3, g3, A, F);   x1 = xx1;   x2 = xx2;   x3 = xx3;}static void StripZeroes(vec_ZZ_p& x){   long n = x.length();   while (n > 0 && IsZero(x[n-1]))      n--;   x.SetLength(n);}void PlainUpdateMap(vec_ZZ_p& xx, const vec_ZZ_p& a,                     const ZZ_pX& b, const ZZ_pX& f){   long n = deg(f);   long i, m;   if (IsZero(b)) {      xx.SetLength(0);      return;   }   m = n-1 - deg(b);   vec_ZZ_p x(INIT_SIZE, n);   for (i = 0; i <= m; i++)      InnerProduct(x[i], a, b.rep, i);   if (deg(b) != 0) {      ZZ_pX c(INIT_SIZE, n);      LeftShift(c, b, m);      for (i = m+1; i < n; i++) {         MulByXMod(c, c, f);         InnerProduct(x[i], a, c.rep);      }   }   xx = x;}   void UpdateMap(vec_ZZ_p& x, const vec_ZZ_p& aa,                const ZZ_pXMultiplier& B, const ZZ_pXModulus& F){   long n = F.n;   long i;   vec_ZZ_p a;   a = aa;   StripZeroes(a);   if (a.length() > n) Error("UpdateMap: bad args");   if (!B.UseFFT) {      PlainUpdateMap(x, a, B.b, F.f);      StripZeroes(x);      return;   }   FFTRep R1(INIT_SIZE, F.k), R2(INIT_SIZE, F.l);   vec_ZZ_p V1(INIT_SIZE, n);   RevToFFTRep(R1, a, F.k, 0, a.length()-1, 0);   mul(R2, R1, F.FRep);   RevFromFFTRep(V1, R2, 0, n-2);   for (i = 0; i <= n-2; i++)  negate(V1[i], V1[i]);   RevToFFTRep(R2, V1, F.l, 0, n-2, n-1);   mul(R2, R2, B.B1);   mul(R1, R1, B.B2);   AddExpand(R2, R1);   RevFromFFTRep(x, R2, 0, n-1);   StripZeroes(x);}   void ProjectPowers(vec_ZZ_p& x, const vec_ZZ_p& a, long k,                   const ZZ_pXArgument& H, const ZZ_pXModulus& F){   long n = F.n;   if (a.length() > n || k < 0 || k >= (1L << (NTL_BITS_PER_LONG-4)))       Error("ProjectPowers: bad args");   long m = H.H.length()-1;   long l = (k+m-1)/m - 1;   ZZ_pXMultiplier M;   build(M, H.H[m], F);   vec_ZZ_p s(INIT_SIZE, n);   s = a;   StripZeroes(s);   x.SetLength(k);   for (long i = 0; i <= l; i++) {      long m1 = min(m, k-i*m);      ZZ_p* w = &x[i*m];      for (long j = 0; j < m1; j++)         InnerProduct(w[j], H.H[j].rep, s);      if (i < l)         UpdateMap(s, s, M, F);   }}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 BerlekampMassey(ZZ_pX& h, const vec_ZZ_p& a, long m){   ZZ_pX Lambda, Sigma, Temp;   long L;   ZZ_p Delta, Delta1, t1;   long shamt;   // cerr << "*** " << m << "\n";   Lambda.SetMaxLength(m+1);   Sigma.SetMaxLength(m+1);   Temp.SetMaxLength(m+1);   L = 0;   set(Lambda);   clear(Sigma);   set(Delta);   shamt = 0;   long i, r, dl;   for (r = 1; r <= 2*m; r++) {      // cerr << r << "--";      clear(Delta1);      dl = deg(Lambda);      for (i = 0; i <= dl; i++) {         mul(t1, Lambda.rep[i], a[r-i-1]);         add(Delta1, Delta1, t1);      }      if (IsZero(Delta1)) {         shamt++;         // cerr << "case 1: " << deg(Lambda) << " " << deg(Sigma) << " " << shamt << "\n";      }      else if (2*L < r) {         div(t1, Delta1, Delta);         mul(Temp, Sigma, t1);         Sigma = Lambda;         ShiftSub(Lambda, Temp, shamt+1);         shamt = 0;         L = r-L;         Delta = Delta1;         // cerr << "case 2: " << deg(Lambda) << " " << deg(Sigma) << " " << shamt << "\n";      }      else {         shamt++;         div(t1, Delta1, Delta);         mul(Temp, Sigma, t1);         ShiftSub(Lambda, Temp, shamt);         // cerr << "case 3: " << deg(Lambda) << " " << deg(Sigma) << " " << shamt << "\n";      }   }   // cerr << "finished: " << L << " " << deg(Lambda) << "\n";    dl = deg(Lambda);   h.rep.SetLength(L + 1);   for (i = 0; i < L - dl; i++)      clear(h.rep[i]);   for (i = L - dl; i <= L; i++)      h.rep[i] = Lambda.rep[L - i];}void GCDMinPolySeq(ZZ_pX& h, const vec_ZZ_p& x, long m){   long i;   ZZ_pX a, b;   ZZ_pXMatrix M;   ZZ_p t;   a.rep.SetLength(2*m);   for (i = 0; i < 2*m; i++) a.rep[i] = x[2*m-1-i];   a.normalize();   SetCoeff(b, 2*m);   HalfGCD(M, b, a, m+1);   /* make monic */   inv(t, LeadCoeff(M(1,1)));   mul(h, M(1,1), t);}void MinPolySeq(ZZ_pX& h, const vec_ZZ_p& a, long m){   if (m < 0 || m >= (1L << (NTL_BITS_PER_LONG-4))) 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 DoMinPolyMod(ZZ_pX& h, const ZZ_pX& g, const ZZ_pXModulus& F, long m,               const vec_ZZ_p& R) {   vec_ZZ_p x;   ProjectPowers(x, R, 2*m, g, F);   MinPolySeq(h, x, m);}