source: git/ntl/src/ZZ_pEX.c @ 33a041

spielwiese
Last change on this file since 33a041 was 311902, checked in by Hans Schönemann <hannes@…>, 20 years ago
*hannes: 5.3.2-C++-fix git-svn-id: file:///usr/local/Singular/svn/trunk@7474 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 57.3 KB
Line 
1
2
3
4#include <NTL/ZZ_pEX.h>
5#include <NTL/vec_vec_ZZ_p.h>
6
7#include <NTL/new.h>
8
9NTL_START_IMPL
10
11
12const ZZ_pEX& ZZ_pEX::zero()
13{
14   static ZZ_pEX z;
15   return z;
16}
17
18
19void ZZ_pEX::normalize()
20{
21   long n;
22   const ZZ_pE* p;
23
24   n = rep.length();
25   if (n == 0) return;
26   p = rep.elts() + n;
27   while (n > 0 && IsZero(*--p)) {
28      n--;
29   }
30   rep.SetLength(n);
31}
32
33
34long IsZero(const ZZ_pEX& a)
35{
36   return a.rep.length() == 0;
37}
38
39
40long IsOne(const ZZ_pEX& a)
41{
42    return a.rep.length() == 1 && IsOne(a.rep[0]);
43}
44
45long operator==(const ZZ_pEX& a, long b)
46{
47   if (b == 0)
48      return IsZero(a);
49
50   if (b == 1)
51      return IsOne(a);
52
53   long da = deg(a);
54
55   if (da > 0) return 0;
56
57   NTL_ZZ_pRegister(bb);
58   bb = b;
59
60   if (da < 0)
61      return IsZero(bb);
62
63   return a.rep[0] == bb;
64}
65
66long operator==(const ZZ_pEX& a, const ZZ_p& b)
67{
68   if (IsZero(b))
69      return IsZero(a);
70
71   long da = deg(a);
72
73   if (da != 0)
74      return 0;
75
76   return a.rep[0] == b;
77}
78
79long operator==(const ZZ_pEX& a, const ZZ_pE& b)
80{
81   if (IsZero(b))
82      return IsZero(a);
83
84   long da = deg(a);
85
86   if (da != 0)
87      return 0;
88
89   return a.rep[0] == b;
90}
91
92
93
94
95void SetCoeff(ZZ_pEX& x, long i, const ZZ_pE& a)
96{
97   long j, m;
98
99   if (i < 0) 
100      Error("SetCoeff: negative index");
101
102   if (NTL_OVERFLOW(i, 1, 0))
103      Error("overflow in SetCoeff");
104
105   m = deg(x);
106
107   if (i > m) {
108      /* careful: a may alias a coefficient of x */
109
110      long alloc = x.rep.allocated();
111
112      if (alloc > 0 && i >= alloc) {
113         ZZ_pE aa = a;
114         x.rep.SetLength(i+1);
115         x.rep[i] = aa;
116      }
117      else {
118         x.rep.SetLength(i+1);
119         x.rep[i] = a;
120      }
121
122      for (j = m+1; j < i; j++)
123         clear(x.rep[j]);
124   }
125   else
126      x.rep[i] = a;
127
128   x.normalize();
129}
130
131void SetCoeff(ZZ_pEX& x, long i, const ZZ_p& aa)
132{
133   long j, m;
134
135   if (i < 0)
136      Error("SetCoeff: negative index");
137
138   if (NTL_OVERFLOW(i, 1, 0))
139      Error("overflow in SetCoeff");
140
141   NTL_ZZ_pRegister(a);  // watch out for aliases!
142   a = aa;
143
144   m = deg(x);
145
146   if (i > m) {
147      x.rep.SetLength(i+1);
148      for (j = m+1; j < i; j++)
149         clear(x.rep[j]);
150   }
151   x.rep[i] = a;
152   x.normalize();
153}
154
155void SetCoeff(ZZ_pEX& x, long i, long a)
156{
157   if (a == 1)
158      SetCoeff(x, i);
159   else {
160      NTL_ZZ_pRegister(T);
161      T = a;
162      SetCoeff(x, i, T);
163   }
164}
165
166
167
168void SetCoeff(ZZ_pEX& x, long i)
169{
170   long j, m;
171
172   if (i < 0) 
173      Error("coefficient index out of range");
174
175   if (NTL_OVERFLOW(i, 1, 0))
176      Error("overflow in SetCoeff");
177
178   m = deg(x);
179
180   if (i > m) {
181      x.rep.SetLength(i+1);
182      for (j = m+1; j < i; j++)
183         clear(x.rep[j]);
184   }
185   set(x.rep[i]);
186   x.normalize();
187}
188
189
190void SetX(ZZ_pEX& x)
191{
192   clear(x);
193   SetCoeff(x, 1);
194}
195
196
197long IsX(const ZZ_pEX& a)
198{
199   return deg(a) == 1 && IsOne(LeadCoeff(a)) && IsZero(ConstTerm(a));
200}
201     
202     
203
204const ZZ_pE& coeff(const ZZ_pEX& a, long i)
205{
206   if (i < 0 || i > deg(a))
207      return ZZ_pE::zero();
208   else
209      return a.rep[i];
210}
211
212
213const ZZ_pE& LeadCoeff(const ZZ_pEX& a)
214{
215   if (IsZero(a))
216      return ZZ_pE::zero();
217   else
218      return a.rep[deg(a)];
219}
220
221const ZZ_pE& ConstTerm(const ZZ_pEX& a)
222{
223   if (IsZero(a))
224      return ZZ_pE::zero();
225   else
226      return a.rep[0];
227}
228
229
230
231void conv(ZZ_pEX& x, const ZZ_pE& a)
232{
233   if (IsZero(a))
234      x.rep.SetLength(0);
235   else {
236      x.rep.SetLength(1);
237      x.rep[0] = a;
238   }
239}
240
241void conv(ZZ_pEX& x, long a)
242{
243   if (a == 0) 
244      clear(x);
245   else if (a == 1)
246      set(x);
247   else {
248      NTL_ZZ_pRegister(T);
249      T = a;
250      conv(x, T);
251   }
252}
253
254void conv(ZZ_pEX& x, const ZZ& a)
255{
256   NTL_ZZ_pRegister(T);
257   conv(T, a);
258   conv(x, T);
259}
260
261void conv(ZZ_pEX& x, const ZZ_p& a)
262{
263   if (IsZero(a)) 
264      clear(x);
265   else if (IsOne(a))
266      set(x);
267   else {
268      x.rep.SetLength(1);
269      conv(x.rep[0], a);
270      x.normalize();
271   }
272}
273
274void conv(ZZ_pEX& x, const ZZ_pX& aa)
275{
276   ZZ_pX a = aa; // in case a aliases the rep of a coefficient of x
277
278   long n = deg(a)+1;
279   long i;
280
281   x.rep.SetLength(n);
282   for (i = 0; i < n; i++)
283      conv(x.rep[i], coeff(a, i));
284}
285
286
287void conv(ZZ_pEX& x, const vec_ZZ_pE& a)
288{
289   x.rep = a;
290   x.normalize();
291}
292
293
294void add(ZZ_pEX& x, const ZZ_pEX& a, const ZZ_pEX& b)
295{
296   long da = deg(a);
297   long db = deg(b);
298   long minab = min(da, db);
299   long maxab = max(da, db);
300   x.rep.SetLength(maxab+1);
301
302   long i;
303   const ZZ_pE *ap, *bp; 
304   ZZ_pE* xp;
305
306   for (i = minab+1, ap = a.rep.elts(), bp = b.rep.elts(), xp = x.rep.elts();
307        i; i--, ap++, bp++, xp++)
308      add(*xp, (*ap), (*bp));
309
310   if (da > minab && &x != &a)
311      for (i = da-minab; i; i--, xp++, ap++)
312         *xp = *ap;
313   else if (db > minab && &x != &b)
314      for (i = db-minab; i; i--, xp++, bp++)
315         *xp = *bp;
316   else
317      x.normalize();
318}
319
320
321void add(ZZ_pEX& x, const ZZ_pEX& a, const ZZ_pE& b)
322{
323   long n = a.rep.length();
324   if (n == 0) {
325      conv(x, b);
326   }
327   else if (&x == &a) {
328      add(x.rep[0], a.rep[0], b);
329      x.normalize();
330   }
331   else if (x.rep.MaxLength() == 0) {
332      x = a;
333      add(x.rep[0], a.rep[0], b);
334      x.normalize();
335   }
336   else {
337      // ugly...b could alias a coeff of x
338
339      ZZ_pE *xp = x.rep.elts();
340      add(xp[0], a.rep[0], b);
341      x.rep.SetLength(n);
342      xp = x.rep.elts();
343      const ZZ_pE *ap = a.rep.elts();
344      long i;
345      for (i = 1; i < n; i++)
346         xp[i] = ap[i];
347      x.normalize();
348   }
349}
350
351void add(ZZ_pEX& x, const ZZ_pEX& a, const ZZ_p& b)
352{
353   long n = a.rep.length();
354   if (n == 0) {
355      conv(x, b);
356   }
357   else if (&x == &a) {
358      add(x.rep[0], a.rep[0], b);
359      x.normalize();
360   }
361   else if (x.rep.MaxLength() == 0) {
362      x = a;
363      add(x.rep[0], a.rep[0], b);
364      x.normalize();
365   }
366   else {
367      // ugly...b could alias a coeff of x
368
369      ZZ_pE *xp = x.rep.elts();
370      add(xp[0], a.rep[0], b);
371      x.rep.SetLength(n);
372      xp = x.rep.elts();
373      const ZZ_pE *ap = a.rep.elts();
374      long i;
375      for (i = 1; i < n; i++)
376         xp[i] = ap[i];
377      x.normalize();
378   }
379}
380
381
382void add(ZZ_pEX& x, const ZZ_pEX& a, long b)
383{
384   if (a.rep.length() == 0) {
385      conv(x, b);
386   }
387   else {
388      if (&x != &a) x = a;
389      add(x.rep[0], x.rep[0], b);
390      x.normalize();
391   }
392}
393
394
395void sub(ZZ_pEX& x, const ZZ_pEX& a, const ZZ_pEX& b)
396{
397   long da = deg(a);
398   long db = deg(b);
399   long minab = min(da, db);
400   long maxab = max(da, db);
401   x.rep.SetLength(maxab+1);
402
403   long i;
404   const ZZ_pE *ap, *bp; 
405   ZZ_pE* xp;
406
407   for (i = minab+1, ap = a.rep.elts(), bp = b.rep.elts(), xp = x.rep.elts();
408        i; i--, ap++, bp++, xp++)
409      sub(*xp, (*ap), (*bp));
410
411   if (da > minab && &x != &a)
412      for (i = da-minab; i; i--, xp++, ap++)
413         *xp = *ap;
414   else if (db > minab)
415      for (i = db-minab; i; i--, xp++, bp++)
416         negate(*xp, *bp);
417   else
418      x.normalize();
419}
420
421
422void sub(ZZ_pEX& x, const ZZ_pEX& a, const ZZ_pE& b)
423{
424   long n = a.rep.length();
425   if (n == 0) {
426      conv(x, b);
427      negate(x, x);
428   }
429   else if (&x == &a) {
430      sub(x.rep[0], a.rep[0], b);
431      x.normalize();
432   }
433   else if (x.rep.MaxLength() == 0) {
434      x = a;
435      sub(x.rep[0], a.rep[0], b);
436      x.normalize();
437   }
438   else {
439      // ugly...b could alias a coeff of x
440
441      ZZ_pE *xp = x.rep.elts();
442      sub(xp[0], a.rep[0], b);
443      x.rep.SetLength(n);
444      xp = x.rep.elts();
445      const ZZ_pE *ap = a.rep.elts();
446      long i;
447      for (i = 1; i < n; i++)
448         xp[i] = ap[i];
449      x.normalize();
450   }
451}
452
453void sub(ZZ_pEX& x, const ZZ_pEX& a, const ZZ_p& b)
454{
455   long n = a.rep.length();
456   if (n == 0) {
457      conv(x, b);
458      negate(x, x);
459   }
460   else if (&x == &a) {
461      sub(x.rep[0], a.rep[0], b);
462      x.normalize();
463   }
464   else if (x.rep.MaxLength() == 0) {
465      x = a;
466      sub(x.rep[0], a.rep[0], b);
467      x.normalize();
468   }
469   else {
470      // ugly...b could alias a coeff of x
471
472      ZZ_pE *xp = x.rep.elts();
473      sub(xp[0], a.rep[0], b);
474      x.rep.SetLength(n);
475      xp = x.rep.elts();
476      const ZZ_pE *ap = a.rep.elts();
477      long i;
478      for (i = 1; i < n; i++)
479         xp[i] = ap[i];
480      x.normalize();
481   }
482}
483
484
485void sub(ZZ_pEX& x, const ZZ_pEX& a, long b)
486{
487   if (a.rep.length() == 0) {
488      conv(x, b);
489      negate(x, x);
490   }
491   else {
492      if (&x != &a) x = a;
493      sub(x.rep[0], x.rep[0], b);
494      x.normalize();
495   }
496}
497
498void sub(ZZ_pEX& x, const ZZ_pE& b, const ZZ_pEX& a)
499{
500   long n = a.rep.length();
501   if (n == 0) {
502      conv(x, b);
503   }
504   else if (x.rep.MaxLength() == 0) {
505      negate(x, a);
506      add(x.rep[0], a.rep[0], b);
507      x.normalize();
508   }
509   else {
510      // ugly...b could alias a coeff of x
511
512      ZZ_pE *xp = x.rep.elts();
513      sub(xp[0], b, a.rep[0]);
514      x.rep.SetLength(n);
515      xp = x.rep.elts();
516      const ZZ_pE *ap = a.rep.elts();
517      long i;
518      for (i = 1; i < n; i++)
519         negate(xp[i], ap[i]);
520      x.normalize();
521   }
522}
523
524
525void sub(ZZ_pEX& x, const ZZ_p& a, const ZZ_pEX& b)
526{
527   NTL_ZZ_pRegister(T);   // avoids aliasing problems
528   T = a;
529   negate(x, b);
530   add(x, x, T);
531}
532
533void sub(ZZ_pEX& x, long a, const ZZ_pEX& b)
534{
535   NTL_ZZ_pRegister(T); 
536   T = a;
537   negate(x, b);
538   add(x, x, T);
539}
540
541void mul(ZZ_pEX& c, const ZZ_pEX& a, const ZZ_pEX& b)
542{
543   if (&a == &b) {
544      sqr(c, a);
545      return;
546   }
547
548   if (IsZero(a) || IsZero(b)) {
549      clear(c);
550      return;
551   }
552
553   if (deg(a) == 0) {
554      mul(c, b, ConstTerm(a));
555      return;
556   } 
557
558   if (deg(b) == 0) {
559      mul(c, a, ConstTerm(b));
560      return;
561   }
562
563   // general case...Kronecker subst
564
565   ZZ_pX A, B, C;
566
567   long da = deg(a);
568   long db = deg(b);
569
570   long n = ZZ_pE::degree();
571   long n2 = 2*n-1;
572
573   if (NTL_OVERFLOW(da+db+1, n2, 0))
574      Error("overflow in ZZ_pEX mul");
575
576   long i, j;
577
578   A.rep.SetLength((da+1)*n2);
579
580   for (i = 0; i <= da; i++) {
581      const ZZ_pX& coeff = rep(a.rep[i]);
582      long dcoeff = deg(coeff);
583      for (j = 0; j <= dcoeff; j++)
584         A.rep[n2*i + j] = coeff.rep[j]; 
585   }
586
587   A.normalize();
588
589   B.rep.SetLength((db+1)*n2);
590
591   for (i = 0; i <= db; i++) {
592      const ZZ_pX& coeff = rep(b.rep[i]);
593      long dcoeff = deg(coeff);
594      for (j = 0; j <= dcoeff; j++)
595         B.rep[n2*i + j] = coeff.rep[j]; 
596   }
597
598   B.normalize();
599
600   mul(C, A, B);
601
602   long Clen = C.rep.length();
603   long lc = (Clen + n2 - 1)/n2;
604   long dc = lc - 1;
605
606   c.rep.SetLength(dc+1);
607
608   ZZ_pX tmp;
609   
610   for (i = 0; i <= dc; i++) {
611      tmp.rep.SetLength(n2);
612      for (j = 0; j < n2 && n2*i + j < Clen; j++)
613         tmp.rep[j] = C.rep[n2*i + j];
614      for (; j < n2; j++)
615         clear(tmp.rep[j]);
616      tmp.normalize();
617      conv(c.rep[i], tmp);
618   }
619 
620   c.normalize();
621}
622
623
624void mul(ZZ_pEX& x, const ZZ_pEX& a, const ZZ_pE& b)
625{
626   if (IsZero(b)) {
627      clear(x);
628      return;
629   }
630
631   ZZ_pE t;
632   t = b;
633
634   long i, da;
635
636   const ZZ_pE *ap;
637   ZZ_pE* xp;
638
639   da = deg(a);
640   x.rep.SetLength(da+1);
641   ap = a.rep.elts();
642   xp = x.rep.elts();
643
644   for (i = 0; i <= da; i++)
645      mul(xp[i], ap[i], t);
646
647   x.normalize();
648}
649
650
651
652void mul(ZZ_pEX& x, const ZZ_pEX& a, const ZZ_p& b)
653{
654   if (IsZero(b)) {
655      clear(x);
656      return;
657   }
658
659   NTL_ZZ_pRegister(t);
660   t = b;
661
662   long i, da;
663
664   const ZZ_pE *ap;
665   ZZ_pE* xp;
666
667   da = deg(a);
668   x.rep.SetLength(da+1);
669   ap = a.rep.elts();
670   xp = x.rep.elts();
671
672   for (i = 0; i <= da; i++)
673      mul(xp[i], ap[i], t);
674
675   x.normalize();
676}
677
678
679void mul(ZZ_pEX& x, const ZZ_pEX& a, long b)
680{
681   NTL_ZZ_pRegister(t);
682   t = b;
683   mul(x, a, t);
684}
685
686void sqr(ZZ_pEX& c, const ZZ_pEX& a)
687{
688   if (IsZero(a)) {
689      clear(c);
690      return;
691   }
692
693   if (deg(a) == 0) {
694      ZZ_pE res;
695      sqr(res, ConstTerm(a));
696      conv(c, res);
697      return;
698   } 
699
700   // general case...Kronecker subst
701
702   ZZ_pX A, C;
703
704   long da = deg(a);
705
706   long n = ZZ_pE::degree();
707   long n2 = 2*n-1;
708
709   if (NTL_OVERFLOW(2*da+1, n2, 0))
710      Error("overflow in ZZ_pEX sqr");
711
712   long i, j;
713
714   A.rep.SetLength((da+1)*n2);
715
716   for (i = 0; i <= da; i++) {
717      const ZZ_pX& coeff = rep(a.rep[i]);
718      long dcoeff = deg(coeff);
719      for (j = 0; j <= dcoeff; j++)
720         A.rep[n2*i + j] = coeff.rep[j]; 
721   }
722
723   A.normalize();
724
725   sqr(C, A);
726
727   long Clen = C.rep.length();
728   long lc = (Clen + n2 - 1)/n2;
729   long dc = lc - 1;
730
731   c.rep.SetLength(dc+1);
732
733   ZZ_pX tmp;
734   
735   for (i = 0; i <= dc; i++) {
736      tmp.rep.SetLength(n2);
737      for (j = 0; j < n2 && n2*i + j < Clen; j++)
738         tmp.rep[j] = C.rep[n2*i + j];
739      for (; j < n2; j++)
740         clear(tmp.rep[j]);
741      tmp.normalize();
742      conv(c.rep[i], tmp);
743   }
744 
745 
746   c.normalize();
747}
748
749
750void MulTrunc(ZZ_pEX& x, const ZZ_pEX& a, const ZZ_pEX& b, long n)
751{
752   if (n < 0) Error("MulTrunc: bad args");
753
754   ZZ_pEX t;
755   mul(t, a, b);
756   trunc(x, t, n);
757}
758
759void SqrTrunc(ZZ_pEX& x, const ZZ_pEX& a, long n)
760{
761   if (n < 0) Error("SqrTrunc: bad args");
762
763   ZZ_pEX t;
764   sqr(t, a);
765   trunc(x, t, n);
766}
767
768
769void CopyReverse(ZZ_pEX& x, const ZZ_pEX& a, long hi)
770
771   // x[0..hi] = reverse(a[0..hi]), with zero fill
772   // input may not alias output
773
774{
775   long i, j, n, m;
776
777   n = hi+1;
778   m = a.rep.length();
779
780   x.rep.SetLength(n);
781
782   const ZZ_pE* ap = a.rep.elts();
783   ZZ_pE* xp = x.rep.elts();
784
785   for (i = 0; i < n; i++) {
786      j = hi-i;
787      if (j < 0 || j >= m)
788         clear(xp[i]);
789      else
790         xp[i] = ap[j];
791   }
792
793   x.normalize();
794} 
795
796
797void trunc(ZZ_pEX& x, const ZZ_pEX& a, long m)
798
799// x = a % X^m, output may alias input
800
801{
802   if (m < 0) Error("trunc: bad args");
803
804   if (&x == &a) {
805      if (x.rep.length() > m) {
806         x.rep.SetLength(m);
807         x.normalize();
808      }
809   }
810   else {
811      long n;
812      long i;
813      ZZ_pE* xp;
814      const ZZ_pE* ap;
815
816      n = min(a.rep.length(), m);
817      x.rep.SetLength(n);
818
819      xp = x.rep.elts();
820      ap = a.rep.elts();
821
822      for (i = 0; i < n; i++) xp[i] = ap[i];
823
824      x.normalize();
825   }
826}
827
828
829void random(ZZ_pEX& x, long n)
830{
831   long i;
832
833   x.rep.SetLength(n);
834
835   for (i = 0; i < n; i++)
836      random(x.rep[i]);
837
838   x.normalize();
839}
840
841void negate(ZZ_pEX& x, const ZZ_pEX& a)
842{
843   long n = a.rep.length();
844   x.rep.SetLength(n);
845
846   const ZZ_pE* ap = a.rep.elts();
847   ZZ_pE* xp = x.rep.elts();
848   long i;
849
850   for (i = n; i; i--, ap++, xp++)
851      negate((*xp), (*ap));
852}
853
854
855
856static
857void MulByXModAux(ZZ_pEX& h, const ZZ_pEX& a, const ZZ_pEX& f)
858{
859   long i, n, m;
860   ZZ_pE* hh;
861   const ZZ_pE *aa, *ff;
862
863   ZZ_pE t, z;
864
865   n = deg(f);
866   m = deg(a);
867
868   if (m >= n || n == 0) Error("MulByXMod: bad args");
869
870   if (m < 0) {
871      clear(h);
872      return;
873   }
874
875   if (m < n-1) {
876      h.rep.SetLength(m+2);
877      hh = h.rep.elts();
878      aa = a.rep.elts();
879      for (i = m+1; i >= 1; i--)
880         hh[i] = aa[i-1];
881      clear(hh[0]);
882   }
883   else {
884      h.rep.SetLength(n);
885      hh = h.rep.elts();
886      aa = a.rep.elts();
887      ff = f.rep.elts();
888      negate(z, aa[n-1]);
889      if (!IsOne(ff[n]))
890         div(z, z, ff[n]);
891      for (i = n-1; i >= 1; i--) {
892         mul(t, z, ff[i]);
893         add(hh[i], aa[i-1], t);
894      }
895      mul(hh[0], z, ff[0]);
896      h.normalize();
897   }
898}
899
900void MulByXMod(ZZ_pEX& h, const ZZ_pEX& a, const ZZ_pEX& f)
901{
902   if (&h == &f) {
903      ZZ_pEX hh;
904      MulByXModAux(hh, a, f);
905      h = hh;
906   }
907   else
908      MulByXModAux(h, a, f);
909}
910
911
912
913void PlainMul(ZZ_pEX& x, const ZZ_pEX& a, const ZZ_pEX& b)
914{
915   long da = deg(a);
916   long db = deg(b);
917
918   if (da < 0 || db < 0) {
919      clear(x);
920      return;
921   }
922
923   long d = da+db;
924
925
926
927   const ZZ_pE *ap, *bp;
928   ZZ_pE *xp;
929   
930   ZZ_pEX la, lb;
931
932   if (&x == &a) {
933      la = a;
934      ap = la.rep.elts();
935   }
936   else
937      ap = a.rep.elts();
938
939   if (&x == &b) {
940      lb = b;
941      bp = lb.rep.elts();
942   }
943   else
944      bp = b.rep.elts();
945
946   x.rep.SetLength(d+1);
947
948   xp = x.rep.elts();
949
950   long i, j, jmin, jmax;
951   static ZZ_pX t, accum;
952
953   for (i = 0; i <= d; i++) {
954      jmin = max(0, i-db);
955      jmax = min(da, i);
956      clear(accum);
957      for (j = jmin; j <= jmax; j++) {
958         mul(t, rep(ap[j]), rep(bp[i-j]));
959         add(accum, accum, t);
960      }
961      conv(xp[i], accum);
962   }
963   x.normalize();
964}
965
966void SetSize(vec_ZZ_pX& x, long n, long m)
967{
968   x.SetLength(n);
969   long i;
970   for (i = 0; i < n; i++)
971      x[i].rep.SetMaxLength(m);
972}
973
974
975
976void PlainDivRem(ZZ_pEX& q, ZZ_pEX& r, const ZZ_pEX& a, const ZZ_pEX& b)
977{
978   long da, db, dq, i, j, LCIsOne;
979   const ZZ_pE *bp;
980   ZZ_pE *qp;
981   ZZ_pX *xp;
982
983
984   ZZ_pE LCInv, t;
985   ZZ_pX s;
986
987   da = deg(a);
988   db = deg(b);
989
990   if (db < 0) Error("ZZ_pEX: division by zero");
991
992   if (da < db) {
993      r = a;
994      clear(q);
995      return;
996   }
997
998   ZZ_pEX lb;
999
1000   if (&q == &b) {
1001      lb = b;
1002      bp = lb.rep.elts();
1003   }
1004   else
1005      bp = b.rep.elts();
1006
1007   if (IsOne(bp[db]))
1008      LCIsOne = 1;
1009   else {
1010      LCIsOne = 0;
1011      inv(LCInv, bp[db]);
1012   }
1013
1014   vec_ZZ_pX x;
1015
1016   SetSize(x, da+1, 2*ZZ_pE::degree());
1017
1018   for (i = 0; i <= da; i++) 
1019      x[i] = rep(a.rep[i]);
1020
1021   xp = x.elts();
1022
1023   dq = da - db;
1024   q.rep.SetLength(dq+1);
1025   qp = q.rep.elts();
1026
1027   for (i = dq; i >= 0; i--) {
1028      conv(t, xp[i+db]);
1029      if (!LCIsOne)
1030         mul(t, t, LCInv);
1031      qp[i] = t;
1032      negate(t, t);
1033
1034      for (j = db-1; j >= 0; j--) {
1035         mul(s, rep(t), rep(bp[j]));
1036         add(xp[i+j], xp[i+j], s);
1037      }
1038   }
1039
1040   r.rep.SetLength(db);
1041   for (i = 0; i < db; i++)
1042      conv(r.rep[i], xp[i]);
1043   r.normalize();
1044}
1045
1046
1047void PlainRem(ZZ_pEX& r, const ZZ_pEX& a, const ZZ_pEX& b, vec_ZZ_pX& x)
1048{
1049   long da, db, dq, i, j, LCIsOne;
1050   const ZZ_pE *bp;
1051   ZZ_pX *xp;
1052
1053
1054   ZZ_pE LCInv, t;
1055   ZZ_pX s;
1056
1057   da = deg(a);
1058   db = deg(b);
1059
1060   if (db < 0) Error("ZZ_pEX: division by zero");
1061
1062   if (da < db) {
1063      r = a;
1064      return;
1065   }
1066
1067   bp = b.rep.elts();
1068
1069   if (IsOne(bp[db]))
1070      LCIsOne = 1;
1071   else {
1072      LCIsOne = 0;
1073      inv(LCInv, bp[db]);
1074   }
1075
1076   for (i = 0; i <= da; i++)
1077      x[i] = rep(a.rep[i]);
1078
1079   xp = x.elts();
1080
1081   dq = da - db;
1082
1083   for (i = dq; i >= 0; i--) {
1084      conv(t, xp[i+db]);
1085      if (!LCIsOne)
1086         mul(t, t, LCInv);
1087      negate(t, t);
1088
1089      for (j = db-1; j >= 0; j--) {
1090         mul(s, rep(t), rep(bp[j]));
1091         add(xp[i+j], xp[i+j], s);
1092      }
1093   }
1094
1095   r.rep.SetLength(db);
1096   for (i = 0; i < db; i++)
1097      conv(r.rep[i], xp[i]);
1098   r.normalize();
1099}
1100
1101
1102void PlainDivRem(ZZ_pEX& q, ZZ_pEX& r, const ZZ_pEX& a, const ZZ_pEX& b, 
1103     vec_ZZ_pX& x)
1104{
1105   long da, db, dq, i, j, LCIsOne;
1106   const ZZ_pE *bp;
1107   ZZ_pE *qp;
1108   ZZ_pX *xp;
1109
1110
1111   ZZ_pE LCInv, t;
1112   ZZ_pX s;
1113
1114   da = deg(a);
1115   db = deg(b);
1116
1117   if (db < 0) Error("ZZ_pEX: division by zero");
1118
1119   if (da < db) {
1120      r = a;
1121      clear(q);
1122      return;
1123   }
1124
1125   ZZ_pEX lb;
1126
1127   if (&q == &b) {
1128      lb = b;
1129      bp = lb.rep.elts();
1130   }
1131   else
1132      bp = b.rep.elts();
1133
1134   if (IsOne(bp[db]))
1135      LCIsOne = 1;
1136   else {
1137      LCIsOne = 0;
1138      inv(LCInv, bp[db]);
1139   }
1140
1141   for (i = 0; i <= da; i++)
1142      x[i] = rep(a.rep[i]);
1143
1144   xp = x.elts();
1145
1146   dq = da - db;
1147   q.rep.SetLength(dq+1);
1148   qp = q.rep.elts();
1149
1150   for (i = dq; i >= 0; i--) {
1151      conv(t, xp[i+db]);
1152      if (!LCIsOne)
1153         mul(t, t, LCInv);
1154      qp[i] = t;
1155      negate(t, t);
1156
1157      for (j = db-1; j >= 0; j--) {
1158         mul(s, rep(t), rep(bp[j]));
1159         add(xp[i+j], xp[i+j], s);
1160      }
1161   }
1162
1163   r.rep.SetLength(db);
1164   for (i = 0; i < db; i++)
1165      conv(r.rep[i], xp[i]);
1166   r.normalize();
1167}
1168
1169
1170void PlainDiv(ZZ_pEX& q, const ZZ_pEX& a, const ZZ_pEX& b)
1171{
1172   long da, db, dq, i, j, LCIsOne;
1173   const ZZ_pE *bp;
1174   ZZ_pE *qp;
1175   ZZ_pX *xp;
1176
1177
1178   ZZ_pE LCInv, t;
1179   ZZ_pX s;
1180
1181   da = deg(a);
1182   db = deg(b);
1183
1184   if (db < 0) Error("ZZ_pEX: division by zero");
1185
1186   if (da < db) {
1187      clear(q);
1188      return;
1189   }
1190
1191   ZZ_pEX lb;
1192
1193   if (&q == &b) {
1194      lb = b;
1195      bp = lb.rep.elts();
1196   }
1197   else
1198      bp = b.rep.elts();
1199
1200   if (IsOne(bp[db]))
1201      LCIsOne = 1;
1202   else {
1203      LCIsOne = 0;
1204      inv(LCInv, bp[db]);
1205   }
1206
1207   vec_ZZ_pX x;
1208   SetSize(x, da+1-db, 2*ZZ_pE::degree());
1209
1210   for (i = db; i <= da; i++)
1211      x[i-db] = rep(a.rep[i]);
1212
1213   xp = x.elts();
1214
1215   dq = da - db;
1216   q.rep.SetLength(dq+1);
1217   qp = q.rep.elts();
1218
1219   for (i = dq; i >= 0; i--) {
1220      conv(t, xp[i]);
1221      if (!LCIsOne)
1222         mul(t, t, LCInv);
1223      qp[i] = t;
1224      negate(t, t);
1225
1226      long lastj = max(0, db-i);
1227
1228      for (j = db-1; j >= lastj; j--) {
1229         mul(s, rep(t), rep(bp[j]));
1230         add(xp[i+j-db], xp[i+j-db], s);
1231      }
1232   }
1233}
1234
1235void PlainRem(ZZ_pEX& r, const ZZ_pEX& a, const ZZ_pEX& b)
1236{
1237   long da, db, dq, i, j, LCIsOne;
1238   const ZZ_pE *bp;
1239   ZZ_pX *xp;
1240
1241
1242   ZZ_pE LCInv, t;
1243   ZZ_pX s;
1244
1245   da = deg(a);
1246   db = deg(b);
1247
1248   if (db < 0) Error("ZZ_pEX: division by zero");
1249
1250   if (da < db) {
1251      r = a;
1252      return;
1253   }
1254
1255   bp = b.rep.elts();
1256
1257   if (IsOne(bp[db]))
1258      LCIsOne = 1;
1259   else {
1260      LCIsOne = 0;
1261      inv(LCInv, bp[db]);
1262   }
1263
1264   vec_ZZ_pX x;
1265   SetSize(x, da + 1, 2*ZZ_pE::degree());
1266
1267   for (i = 0; i <= da; i++)
1268      x[i] = rep(a.rep[i]);
1269
1270   xp = x.elts();
1271
1272   dq = da - db;
1273
1274   for (i = dq; i >= 0; i--) {
1275      conv(t, xp[i+db]);
1276      if (!LCIsOne)
1277         mul(t, t, LCInv);
1278      negate(t, t);
1279
1280      for (j = db-1; j >= 0; j--) {
1281         mul(s, rep(t), rep(bp[j]));
1282         add(xp[i+j], xp[i+j], s);
1283      }
1284   }
1285
1286   r.rep.SetLength(db);
1287   for (i = 0; i < db; i++)
1288      conv(r.rep[i], xp[i]);
1289   r.normalize();
1290}
1291
1292
1293
1294void RightShift(ZZ_pEX& x, const ZZ_pEX& a, long n)
1295{
1296   if (IsZero(a)) {
1297      clear(x);
1298      return;
1299   }
1300
1301   if (n < 0) {
1302      if (n < -NTL_MAX_LONG) Error("overflow in RightShift");
1303      LeftShift(x, a, -n);
1304      return;
1305   }
1306
1307   long da = deg(a);
1308   long i;
1309 
1310   if (da < n) {
1311      clear(x);
1312      return;
1313   }
1314
1315   if (&x != &a)
1316      x.rep.SetLength(da-n+1);
1317
1318   for (i = 0; i <= da-n; i++)
1319      x.rep[i] = a.rep[i+n];
1320
1321   if (&x == &a)
1322      x.rep.SetLength(da-n+1);
1323
1324   x.normalize();
1325}
1326
1327void LeftShift(ZZ_pEX& x, const ZZ_pEX& a, long n)
1328{
1329   if (IsZero(a)) {
1330      clear(x);
1331      return;
1332   }
1333
1334   if (n < 0) {
1335      if (n < -NTL_MAX_LONG) 
1336         clear(x);
1337      else
1338         RightShift(x, a, -n);
1339      return;
1340   }
1341
1342   if (NTL_OVERFLOW(n, 1, 0))
1343      Error("overflow in LeftShift");
1344
1345   long m = a.rep.length();
1346
1347   x.rep.SetLength(m+n);
1348
1349   long i;
1350   for (i = m-1; i >= 0; i--)
1351      x.rep[i+n] = a.rep[i];
1352
1353   for (i = 0; i < n; i++)
1354      clear(x.rep[i]);
1355}
1356
1357
1358
1359void NewtonInv(ZZ_pEX& c, const ZZ_pEX& a, long e)
1360{
1361   ZZ_pE x;
1362
1363   inv(x, ConstTerm(a));
1364
1365   if (e == 1) {
1366      conv(c, x);
1367      return;
1368   }
1369
1370   static vec_long E;
1371   E.SetLength(0);
1372   append(E, e);
1373   while (e > 1) {
1374      e = (e+1)/2;
1375      append(E, e);
1376   }
1377
1378   long L = E.length();
1379
1380   ZZ_pEX g, g0, g1, g2;
1381
1382
1383   g.rep.SetMaxLength(E[0]);
1384   g0.rep.SetMaxLength(E[0]);
1385   g1.rep.SetMaxLength((3*E[0]+1)/2);
1386   g2.rep.SetMaxLength(E[0]);
1387
1388   conv(g, x);
1389
1390   long i;
1391
1392   for (i = L-1; i > 0; i--) {
1393      // lift from E[i] to E[i-1]
1394
1395      long k = E[i];
1396      long l = E[i-1]-E[i];
1397
1398      trunc(g0, a, k+l);
1399
1400      mul(g1, g0, g);
1401      RightShift(g1, g1, k);
1402      trunc(g1, g1, l);
1403
1404      mul(g2, g1, g);
1405      trunc(g2, g2, l);
1406      LeftShift(g2, g2, k);
1407
1408      sub(g, g, g2);
1409   }
1410
1411   c = g;
1412}
1413
1414void InvTrunc(ZZ_pEX& c, const ZZ_pEX& a, long e)
1415{
1416   if (e < 0) Error("InvTrunc: bad args");
1417   if (e == 0) {
1418      clear(c);
1419      return;
1420   }
1421
1422   if (NTL_OVERFLOW(e, 1, 0))
1423      Error("overflow in InvTrunc");
1424
1425   NewtonInv(c, a, e);
1426}
1427
1428
1429
1430
1431const long ZZ_pEX_MOD_PLAIN = 0;
1432const long ZZ_pEX_MOD_MUL = 1;
1433
1434
1435void build(ZZ_pEXModulus& F, const ZZ_pEX& f)
1436{
1437   long n = deg(f);
1438
1439   if (n <= 0) Error("build(ZZ_pEXModulus,ZZ_pEX): deg(f) <= 0");
1440
1441   if (NTL_OVERFLOW(n, ZZ_pE::degree(), 0))
1442      Error("build(ZZ_pEXModulus,ZZ_pEX): overflow");
1443   
1444
1445   F.tracevec.SetLength(0);
1446
1447   F.f = f;
1448   F.n = n;
1449
1450   if (F.n < ZZ_pE::ModCross()) {
1451      F.method = ZZ_pEX_MOD_PLAIN;
1452   }
1453   else {
1454      F.method = ZZ_pEX_MOD_MUL;
1455      ZZ_pEX P1;
1456      ZZ_pEX P2;
1457
1458      CopyReverse(P1, f, n);
1459      InvTrunc(P2, P1, n-1);
1460      CopyReverse(P1, P2, n-2);
1461      trunc(F.h0, P1, n-2);
1462      trunc(F.f0, f, n);
1463      F.hlc = ConstTerm(P2);
1464   }
1465}
1466
1467
1468
1469ZZ_pEXModulus::ZZ_pEXModulus()
1470{
1471   n = -1;
1472   method = ZZ_pEX_MOD_PLAIN;
1473}
1474
1475
1476ZZ_pEXModulus::~ZZ_pEXModulus() 
1477{ 
1478}
1479
1480
1481
1482ZZ_pEXModulus::ZZ_pEXModulus(const ZZ_pEX& ff)
1483{
1484   n = -1;
1485   method = ZZ_pEX_MOD_PLAIN;
1486
1487   build(*this, ff);
1488}
1489
1490
1491void UseMulRem21(ZZ_pEX& r, const ZZ_pEX& a, const ZZ_pEXModulus& F)
1492{
1493   ZZ_pEX P1;
1494   ZZ_pEX P2;
1495
1496   RightShift(P1, a, F.n);
1497   mul(P2, P1, F.h0);
1498   RightShift(P2, P2, F.n-2);
1499   if (!IsOne(F.hlc)) mul(P1, P1, F.hlc);
1500   add(P2, P2, P1);
1501   mul(P1, P2, F.f0);
1502   trunc(P1, P1, F.n);
1503   trunc(r, a, F.n);
1504   sub(r, r, P1);
1505}
1506
1507void UseMulDivRem21(ZZ_pEX& q, ZZ_pEX& r, const ZZ_pEX& a, const ZZ_pEXModulus& F)
1508{
1509   ZZ_pEX P1;
1510   ZZ_pEX P2;
1511
1512   RightShift(P1, a, F.n);
1513   mul(P2, P1, F.h0);
1514   RightShift(P2, P2, F.n-2);
1515   if (!IsOne(F.hlc)) mul(P1, P1, F.hlc);
1516   add(P2, P2, P1);
1517   mul(P1, P2, F.f0);
1518   trunc(P1, P1, F.n);
1519   trunc(r, a, F.n);
1520   sub(r, r, P1);
1521   q = P2;
1522}
1523
1524void UseMulDiv21(ZZ_pEX& q, const ZZ_pEX& a, const ZZ_pEXModulus& F)
1525{
1526   ZZ_pEX P1;
1527   ZZ_pEX P2;
1528
1529   RightShift(P1, a, F.n);
1530   mul(P2, P1, F.h0);
1531   RightShift(P2, P2, F.n-2);
1532   if (!IsOne(F.hlc)) mul(P1, P1, F.hlc);
1533   add(P2, P2, P1);
1534   q = P2;
1535
1536}
1537
1538
1539void rem(ZZ_pEX& x, const ZZ_pEX& a, const ZZ_pEXModulus& F)
1540{
1541   if (F.method == ZZ_pEX_MOD_PLAIN) {
1542      PlainRem(x, a, F.f);
1543      return;
1544   }
1545
1546   long da = deg(a);
1547   long n = F.n;
1548
1549   if (da <= 2*n-2) {
1550      UseMulRem21(x, a, F);
1551      return;
1552   }
1553
1554   ZZ_pEX buf(INIT_SIZE, 2*n-1);
1555
1556   long a_len = da+1;
1557
1558   while (a_len > 0) {
1559      long old_buf_len = buf.rep.length();
1560      long amt = min(2*n-1-old_buf_len, a_len);
1561
1562      buf.rep.SetLength(old_buf_len+amt);
1563
1564      long i;
1565
1566      for (i = old_buf_len+amt-1; i >= amt; i--)
1567         buf.rep[i] = buf.rep[i-amt];
1568
1569      for (i = amt-1; i >= 0; i--)
1570         buf.rep[i] = a.rep[a_len-amt+i];
1571
1572      buf.normalize();
1573
1574      UseMulRem21(buf, buf, F);
1575
1576      a_len -= amt;
1577   }
1578
1579   x = buf;
1580}
1581
1582void DivRem(ZZ_pEX& q, ZZ_pEX& r, const ZZ_pEX& a, const ZZ_pEXModulus& F)
1583{
1584   if (F.method == ZZ_pEX_MOD_PLAIN) {
1585      PlainDivRem(q, r, a, F.f);
1586      return;
1587   }
1588
1589   long da = deg(a);
1590   long n = F.n;
1591
1592   if (da <= 2*n-2) {
1593      UseMulDivRem21(q, r, a, F);
1594      return;
1595   }
1596
1597   ZZ_pEX buf(INIT_SIZE, 2*n-1);
1598   ZZ_pEX qbuf(INIT_SIZE, n-1);
1599
1600   ZZ_pEX qq;
1601   qq.rep.SetLength(da-n+1);
1602
1603   long a_len = da+1;
1604   long q_hi = da-n+1;
1605
1606   while (a_len > 0) {
1607      long old_buf_len = buf.rep.length();
1608      long amt = min(2*n-1-old_buf_len, a_len);
1609
1610      buf.rep.SetLength(old_buf_len+amt);
1611
1612      long i;
1613
1614      for (i = old_buf_len+amt-1; i >= amt; i--)
1615         buf.rep[i] = buf.rep[i-amt];
1616
1617      for (i = amt-1; i >= 0; i--)
1618         buf.rep[i] = a.rep[a_len-amt+i];
1619
1620      buf.normalize();
1621
1622      UseMulDivRem21(qbuf, buf, buf, F);
1623      long dl = qbuf.rep.length();
1624      a_len = a_len - amt;
1625      for(i = 0; i < dl; i++)
1626         qq.rep[a_len+i] = qbuf.rep[i];
1627      for(i = dl+a_len; i < q_hi; i++)
1628         clear(qq.rep[i]);
1629      q_hi = a_len;
1630   }
1631
1632   r = buf;
1633
1634   qq.normalize();
1635   q = qq;
1636}
1637
1638void div(ZZ_pEX& q, const ZZ_pEX& a, const ZZ_pEXModulus& F)
1639{
1640   if (F.method == ZZ_pEX_MOD_PLAIN) {
1641      PlainDiv(q, a, F.f);
1642      return;
1643   }
1644
1645   long da = deg(a);
1646   long n = F.n;
1647
1648   if (da <= 2*n-2) {
1649      UseMulDiv21(q, a, F);
1650      return;
1651   }
1652
1653   ZZ_pEX buf(INIT_SIZE, 2*n-1);
1654   ZZ_pEX qbuf(INIT_SIZE, n-1);
1655
1656   ZZ_pEX qq;
1657   qq.rep.SetLength(da-n+1);
1658
1659   long a_len = da+1;
1660   long q_hi = da-n+1;
1661
1662   while (a_len > 0) {
1663      long old_buf_len = buf.rep.length();
1664      long amt = min(2*n-1-old_buf_len, a_len);
1665
1666      buf.rep.SetLength(old_buf_len+amt);
1667
1668      long i;
1669
1670      for (i = old_buf_len+amt-1; i >= amt; i--)
1671         buf.rep[i] = buf.rep[i-amt];
1672
1673      for (i = amt-1; i >= 0; i--)
1674         buf.rep[i] = a.rep[a_len-amt+i];
1675
1676      buf.normalize();
1677
1678      a_len = a_len - amt;
1679      if (a_len > 0)
1680         UseMulDivRem21(qbuf, buf, buf, F);
1681      else
1682         UseMulDiv21(qbuf, buf, F);
1683
1684      long dl = qbuf.rep.length();
1685      for(i = 0; i < dl; i++)
1686         qq.rep[a_len+i] = qbuf.rep[i];
1687      for(i = dl+a_len; i < q_hi; i++)
1688         clear(qq.rep[i]);
1689      q_hi = a_len;
1690   }
1691
1692   qq.normalize();
1693   q = qq;
1694}
1695
1696
1697
1698
1699void MulMod(ZZ_pEX& c, const ZZ_pEX& a, const ZZ_pEX& b, const ZZ_pEXModulus& F)
1700{
1701   if (deg(a) >= F.n || deg(b) >= F.n) Error("MulMod: bad args");
1702
1703   ZZ_pEX t;
1704   mul(t, a, b);
1705   rem(c, t, F);
1706}
1707
1708
1709void SqrMod(ZZ_pEX& c, const ZZ_pEX& a, const ZZ_pEXModulus& F)
1710{
1711   if (deg(a) >= F.n) Error("MulMod: bad args");
1712
1713   ZZ_pEX t;
1714   sqr(t, a);
1715   rem(c, t, F);
1716}
1717
1718
1719
1720void UseMulRem(ZZ_pEX& r, const ZZ_pEX& a, const ZZ_pEX& b)
1721{
1722   ZZ_pEX P1;
1723   ZZ_pEX P2;
1724
1725   long da = deg(a);
1726   long db = deg(b);
1727
1728   CopyReverse(P1, b, db);
1729   InvTrunc(P2, P1, da-db+1);
1730   CopyReverse(P1, P2, da-db);
1731
1732   RightShift(P2, a, db);
1733   mul(P2, P1, P2);
1734   RightShift(P2, P2, da-db);
1735   mul(P1, P2, b);
1736   sub(P1, a, P1);
1737   
1738   r = P1;
1739}
1740
1741void UseMulDivRem(ZZ_pEX& q, ZZ_pEX& r, const ZZ_pEX& a, const ZZ_pEX& b)
1742{
1743   ZZ_pEX P1;
1744   ZZ_pEX P2;
1745
1746   long da = deg(a);
1747   long db = deg(b);
1748
1749   CopyReverse(P1, b, db);
1750   InvTrunc(P2, P1, da-db+1);
1751   CopyReverse(P1, P2, da-db);
1752
1753   RightShift(P2, a, db);
1754   mul(P2, P1, P2);
1755   RightShift(P2, P2, da-db);
1756   mul(P1, P2, b);
1757   sub(P1, a, P1);
1758   
1759   r = P1;
1760   q = P2;
1761}
1762
1763void UseMulDiv(ZZ_pEX& q, const ZZ_pEX& a, const ZZ_pEX& b)
1764{
1765   ZZ_pEX P1;
1766   ZZ_pEX P2;
1767
1768   long da = deg(a);
1769   long db = deg(b);
1770
1771   CopyReverse(P1, b, db);
1772   InvTrunc(P2, P1, da-db+1);
1773   CopyReverse(P1, P2, da-db);
1774
1775   RightShift(P2, a, db);
1776   mul(P2, P1, P2);
1777   RightShift(P2, P2, da-db);
1778   
1779   q = P2;
1780}
1781
1782
1783
1784void DivRem(ZZ_pEX& q, ZZ_pEX& r, const ZZ_pEX& a, const ZZ_pEX& b)
1785{
1786   long sa = a.rep.length();
1787   long sb = b.rep.length();
1788
1789   if (sb < ZZ_pE::DivCross() || sa-sb < ZZ_pE::DivCross())
1790      PlainDivRem(q, r, a, b);
1791   else if (sa < 4*sb)
1792      UseMulDivRem(q, r, a, b);
1793   else {
1794      ZZ_pEXModulus B;
1795      build(B, b);
1796      DivRem(q, r, a, B);
1797   }
1798}
1799
1800void div(ZZ_pEX& q, const ZZ_pEX& a, const ZZ_pEX& b)
1801{
1802   long sa = a.rep.length();
1803   long sb = b.rep.length();
1804
1805   if (sb < ZZ_pE::DivCross() || sa-sb < ZZ_pE::DivCross())
1806      PlainDiv(q, a, b);
1807   else if (sa < 4*sb)
1808      UseMulDiv(q, a, b);
1809   else {
1810      ZZ_pEXModulus B;
1811      build(B, b);
1812      div(q, a, B);
1813   }
1814}
1815
1816void div(ZZ_pEX& q, const ZZ_pEX& a, const ZZ_pE& b)
1817{
1818   ZZ_pE T;
1819   inv(T, b);
1820   mul(q, a, T);
1821}
1822
1823void div(ZZ_pEX& q, const ZZ_pEX& a, const ZZ_p& b)
1824{
1825   NTL_ZZ_pRegister(T);
1826   inv(T, b);
1827   mul(q, a, T);
1828}
1829
1830void div(ZZ_pEX& q, const ZZ_pEX& a, long b)
1831{
1832   NTL_ZZ_pRegister(T);
1833   T = b;
1834   inv(T, T);
1835   mul(q, a, T);
1836}
1837
1838void rem(ZZ_pEX& r, const ZZ_pEX& a, const ZZ_pEX& b)
1839{
1840   long sa = a.rep.length();
1841   long sb = b.rep.length();
1842
1843   if (sb < ZZ_pE::DivCross() || sa-sb < ZZ_pE::DivCross())
1844      PlainRem(r, a, b);
1845   else if (sa < 4*sb)
1846      UseMulRem(r, a, b);
1847   else {
1848      ZZ_pEXModulus B;
1849      build(B, b);
1850      rem(r, a, B);
1851   }
1852}
1853
1854void GCD(ZZ_pEX& x, const ZZ_pEX& a, const ZZ_pEX& b)
1855{
1856   ZZ_pE t;
1857
1858   if (IsZero(b))
1859      x = a;
1860   else if (IsZero(a))
1861      x = b;
1862   else {
1863      long n = max(deg(a),deg(b)) + 1;
1864      ZZ_pEX u(INIT_SIZE, n), v(INIT_SIZE, n);
1865
1866      vec_ZZ_pX tmp;
1867      SetSize(tmp, n, 2*ZZ_pE::degree());
1868
1869      u = a;
1870      v = b;
1871      do {
1872         PlainRem(u, u, v, tmp);
1873         swap(u, v);
1874      } while (!IsZero(v));
1875
1876      x = u;
1877   }
1878
1879   if (IsZero(x)) return;
1880   if (IsOne(LeadCoeff(x))) return;
1881
1882   /* make gcd monic */
1883
1884
1885   inv(t, LeadCoeff(x)); 
1886   mul(x, x, t); 
1887}
1888
1889
1890
1891         
1892
1893void XGCD(ZZ_pEX& d, ZZ_pEX& s, ZZ_pEX& t, const ZZ_pEX& a, const ZZ_pEX& b)
1894{
1895   ZZ_pE z;
1896
1897
1898   if (IsZero(b)) {
1899      set(s);
1900      clear(t);
1901      d = a;
1902   }
1903   else if (IsZero(a)) {
1904      clear(s);
1905      set(t);
1906      d = b;
1907   }
1908   else {
1909      long e = max(deg(a), deg(b)) + 1;
1910
1911      ZZ_pEX temp(INIT_SIZE, e), u(INIT_SIZE, e), v(INIT_SIZE, e), 
1912            u0(INIT_SIZE, e), v0(INIT_SIZE, e), 
1913            u1(INIT_SIZE, e), v1(INIT_SIZE, e), 
1914            u2(INIT_SIZE, e), v2(INIT_SIZE, e), q(INIT_SIZE, e);
1915
1916
1917      set(u1); clear(v1);
1918      clear(u2); set(v2);
1919      u = a; v = b;
1920
1921      do {
1922         DivRem(q, u, u, v);
1923         swap(u, v);
1924         u0 = u2;
1925         v0 = v2;
1926         mul(temp, q, u2);
1927         sub(u2, u1, temp);
1928         mul(temp, q, v2);
1929         sub(v2, v1, temp);
1930         u1 = u0;
1931         v1 = v0;
1932      } while (!IsZero(v));
1933
1934      d = u;
1935      s = u1;
1936      t = v1;
1937   }
1938
1939   if (IsZero(d)) return;
1940   if (IsOne(LeadCoeff(d))) return;
1941
1942   /* make gcd monic */
1943
1944   inv(z, LeadCoeff(d));
1945   mul(d, d, z);
1946   mul(s, s, z);
1947   mul(t, t, z);
1948}
1949
1950NTL_vector_impl(ZZ_pEX,vec_ZZ_pEX)
1951
1952NTL_eq_vector_impl(ZZ_pEX,vec_ZZ_pEX)
1953
1954void IterBuild(ZZ_pE* a, long n)
1955{
1956   long i, k;
1957   ZZ_pE b, t;
1958
1959   if (n <= 0) return;
1960
1961   negate(a[0], a[0]);
1962
1963   for (k = 1; k <= n-1; k++) {
1964      negate(b, a[k]);
1965      add(a[k], b, a[k-1]);
1966      for (i = k-1; i >= 1; i--) {
1967         mul(t, a[i], b);
1968         add(a[i], t, a[i-1]);
1969      }
1970      mul(a[0], a[0], b);
1971   }
1972}
1973
1974void BuildFromRoots(ZZ_pEX& x, const vec_ZZ_pE& a)
1975{
1976   long n = a.length();
1977
1978   if (n == 0) {
1979      set(x);
1980      return;
1981   }
1982
1983   x.rep.SetMaxLength(n+1);
1984   x.rep = a;
1985   IterBuild(&x.rep[0], n);
1986   x.rep.SetLength(n+1);
1987   SetCoeff(x, n);
1988}
1989
1990void eval(ZZ_pE& b, const ZZ_pEX& f, const ZZ_pE& a)
1991// does a Horner evaluation
1992{
1993   ZZ_pE acc;
1994   long i;
1995
1996   clear(acc);
1997   for (i = deg(f); i >= 0; i--) {
1998      mul(acc, acc, a);
1999      add(acc, acc, f.rep[i]);
2000   }
2001
2002   b = acc;
2003}
2004
2005void eval(vec_ZZ_pE& b, const ZZ_pEX& f, const vec_ZZ_pE& a)
2006// naive algorithm:  repeats Horner
2007{
2008   if (&b == &f.rep) {
2009      vec_ZZ_pE bb;
2010      eval(bb, f, a);
2011      b = bb;
2012      return;
2013   }
2014
2015   long m = a.length();
2016   b.SetLength(m);
2017   long i;
2018   for (i = 0; i < m; i++)
2019      eval(b[i], f, a[i]);
2020}
2021
2022
2023void interpolate(ZZ_pEX& f, const vec_ZZ_pE& a, const vec_ZZ_pE& b)
2024{
2025   long m = a.length();
2026   if (b.length() != m) Error("interpolate: vector length mismatch");
2027
2028   if (m == 0) {
2029      clear(f);
2030      return;
2031   }
2032
2033   vec_ZZ_pE prod;
2034   prod = a;
2035
2036   ZZ_pE t1, t2;
2037
2038   long k, i;
2039
2040   vec_ZZ_pE res;
2041   res.SetLength(m);
2042
2043   for (k = 0; k < m; k++) {
2044
2045      const ZZ_pE& aa = a[k];
2046
2047      set(t1);
2048      for (i = k-1; i >= 0; i--) {
2049         mul(t1, t1, aa);
2050         add(t1, t1, prod[i]);
2051      }
2052
2053      clear(t2);
2054      for (i = k-1; i >= 0; i--) {
2055         mul(t2, t2, aa);
2056         add(t2, t2, res[i]);
2057      }
2058
2059
2060      inv(t1, t1);
2061      sub(t2, b[k], t2);
2062      mul(t1, t1, t2);
2063
2064      for (i = 0; i < k; i++) {
2065         mul(t2, prod[i], t1);
2066         add(res[i], res[i], t2);
2067      }
2068
2069      res[k] = t1;
2070
2071      if (k < m-1) {
2072         if (k == 0)
2073            negate(prod[0], prod[0]);
2074         else {
2075            negate(t1, a[k]);
2076            add(prod[k], t1, prod[k-1]);
2077            for (i = k-1; i >= 1; i--) {
2078               mul(t2, prod[i], t1);
2079               add(prod[i], t2, prod[i-1]);
2080            }
2081            mul(prod[0], prod[0], t1);
2082         }
2083      }
2084   }
2085
2086   while (m > 0 && IsZero(res[m-1])) m--;
2087   res.SetLength(m);
2088   f.rep = res;
2089}
2090   
2091void InnerProduct(ZZ_pEX& x, const vec_ZZ_pE& v, long low, long high, 
2092                   const vec_ZZ_pEX& H, long n, vec_ZZ_pX& t)
2093{
2094   ZZ_pX s;
2095   long i, j;
2096
2097   for (j = 0; j < n; j++)
2098      clear(t[j]);
2099
2100   high = min(high, v.length()-1);
2101   for (i = low; i <= high; i++) {
2102      const vec_ZZ_pE& h = H[i-low].rep;
2103      long m = h.length();
2104      const ZZ_pX& w = rep(v[i]);
2105
2106      for (j = 0; j < m; j++) {
2107         mul(s, w, rep(h[j]));
2108         add(t[j], t[j], s);
2109      }
2110   }
2111
2112   x.rep.SetLength(n);
2113   for (j = 0; j < n; j++)
2114      conv(x.rep[j], t[j]);
2115   x.normalize();
2116}
2117
2118
2119
2120void CompMod(ZZ_pEX& x, const ZZ_pEX& g, const ZZ_pEXArgument& A, 
2121             const ZZ_pEXModulus& F)
2122{
2123   if (deg(g) <= 0) {
2124      x = g;
2125      return;
2126   }
2127
2128
2129   ZZ_pEX s, t;
2130   vec_ZZ_pX scratch;
2131   SetSize(scratch, deg(F), 2*ZZ_pE::degree());
2132
2133   long m = A.H.length() - 1;
2134   long l = ((g.rep.length()+m-1)/m) - 1;
2135
2136   const ZZ_pEX& M = A.H[m];
2137
2138   InnerProduct(t, g.rep, l*m, l*m + m - 1, A.H, F.n, scratch);
2139   for (long i = l-1; i >= 0; i--) {
2140      InnerProduct(s, g.rep, i*m, i*m + m - 1, A.H, F.n, scratch);
2141      MulMod(t, t, M, F);
2142      add(t, t, s);
2143   }
2144
2145   x = t;
2146}
2147
2148
2149void build(ZZ_pEXArgument& A, const ZZ_pEX& h, const ZZ_pEXModulus& F, long m)
2150{
2151   long i;
2152
2153   if (m <= 0 || deg(h) >= F.n)
2154      Error("build: bad args");
2155
2156   if (m > F.n) m = F.n;
2157
2158   if (ZZ_pEXArgBound > 0) {
2159      double sz = ZZ_p::storage();
2160      sz = sz*ZZ_pE::degree();
2161      sz = sz + NTL_VECTOR_HEADER_SIZE + sizeof(vec_ZZ_p);
2162      sz = sz*F.n;
2163      sz = sz + NTL_VECTOR_HEADER_SIZE + sizeof(vec_ZZ_pE);
2164      sz = sz/1024;
2165      m = min(m, long(ZZ_pEXArgBound/sz));
2166      m = max(m, 1);
2167   }
2168
2169
2170   A.H.SetLength(m+1);
2171
2172   set(A.H[0]);
2173   A.H[1] = h;
2174   for (i = 2; i <= m; i++)
2175      MulMod(A.H[i], A.H[i-1], h, F);
2176}
2177
2178long ZZ_pEXArgBound = 0;
2179
2180
2181
2182
2183void CompMod(ZZ_pEX& x, const ZZ_pEX& g, const ZZ_pEX& h, const ZZ_pEXModulus& F)
2184   // x = g(h) mod f
2185{
2186   long m = SqrRoot(g.rep.length());
2187
2188   if (m == 0) {
2189      clear(x);
2190      return;
2191   }
2192
2193   ZZ_pEXArgument A;
2194
2195   build(A, h, F, m);
2196
2197   CompMod(x, g, A, F);
2198}
2199
2200
2201
2202
2203void Comp2Mod(ZZ_pEX& x1, ZZ_pEX& x2, const ZZ_pEX& g1, const ZZ_pEX& g2,
2204              const ZZ_pEX& h, const ZZ_pEXModulus& F)
2205
2206{
2207   long m = SqrRoot(g1.rep.length() + g2.rep.length());
2208
2209   if (m == 0) {
2210      clear(x1);
2211      clear(x2);
2212      return;
2213   }
2214
2215   ZZ_pEXArgument A;
2216
2217   build(A, h, F, m);
2218
2219   ZZ_pEX xx1, xx2;
2220
2221   CompMod(xx1, g1, A, F);
2222   CompMod(xx2, g2, A, F);
2223
2224   x1 = xx1;
2225   x2 = xx2;
2226}
2227
2228void Comp3Mod(ZZ_pEX& x1, ZZ_pEX& x2, ZZ_pEX& x3, 
2229              const ZZ_pEX& g1, const ZZ_pEX& g2, const ZZ_pEX& g3,
2230              const ZZ_pEX& h, const ZZ_pEXModulus& F)
2231
2232{
2233   long m = SqrRoot(g1.rep.length() + g2.rep.length() + g3.rep.length());
2234
2235   if (m == 0) {
2236      clear(x1);
2237      clear(x2);
2238      clear(x3);
2239      return;
2240   }
2241
2242   ZZ_pEXArgument A;
2243
2244   build(A, h, F, m);
2245
2246   ZZ_pEX xx1, xx2, xx3;
2247
2248   CompMod(xx1, g1, A, F);
2249   CompMod(xx2, g2, A, F);
2250   CompMod(xx3, g3, A, F);
2251
2252   x1 = xx1;
2253   x2 = xx2;
2254   x3 = xx3;
2255}
2256
2257void build(ZZ_pEXTransMultiplier& B, const ZZ_pEX& b, const ZZ_pEXModulus& F)
2258{
2259   long db = deg(b);
2260
2261   if (db >= F.n) Error("build TransMultiplier: bad args");
2262
2263   ZZ_pEX t;
2264
2265   LeftShift(t, b, F.n-1);
2266   div(t, t, F);
2267
2268   // we optimize for low degree b
2269
2270   long d;
2271
2272   d = deg(t);
2273   if (d < 0)
2274      B.shamt_fbi = 0;
2275   else
2276      B.shamt_fbi = F.n-2 - d; 
2277
2278   CopyReverse(B.fbi, t, d);
2279
2280   // The following code optimizes the case when
2281   // f = X^n + low degree poly
2282
2283   trunc(t, F.f, F.n);
2284   d = deg(t);
2285   if (d < 0)
2286      B.shamt = 0;
2287   else
2288      B.shamt = d;
2289
2290   CopyReverse(B.f0, t, d);
2291
2292   if (db < 0)
2293      B.shamt_b = 0;
2294   else
2295      B.shamt_b = db;
2296
2297   CopyReverse(B.b, b, db);
2298}
2299
2300void TransMulMod(ZZ_pEX& x, const ZZ_pEX& a, const ZZ_pEXTransMultiplier& B,
2301               const ZZ_pEXModulus& F)
2302{
2303   if (deg(a) >= F.n) Error("TransMulMod: bad args");
2304
2305   ZZ_pEX t1, t2;
2306
2307   mul(t1, a, B.b);
2308   RightShift(t1, t1, B.shamt_b);
2309
2310   mul(t2, a, B.f0);
2311   RightShift(t2, t2, B.shamt);
2312   trunc(t2, t2, F.n-1);
2313
2314   mul(t2, t2, B.fbi);
2315   if (B.shamt_fbi > 0) LeftShift(t2, t2, B.shamt_fbi);
2316   trunc(t2, t2, F.n-1);
2317   LeftShift(t2, t2, 1);
2318
2319   sub(x, t1, t2);
2320}
2321
2322
2323void ShiftSub(ZZ_pEX& U, const ZZ_pEX& V, long n)
2324// assumes input does not alias output
2325{
2326   if (IsZero(V))
2327      return;
2328
2329   long du = deg(U);
2330   long dv = deg(V);
2331
2332   long d = max(du, n+dv);
2333
2334   U.rep.SetLength(d+1);
2335   long i;
2336
2337   for (i = du+1; i <= d; i++)
2338      clear(U.rep[i]);
2339
2340   for (i = 0; i <= dv; i++)
2341      sub(U.rep[i+n], U.rep[i+n], V.rep[i]);
2342
2343   U.normalize();
2344}
2345
2346
2347void UpdateMap(vec_ZZ_pE& x, const vec_ZZ_pE& a,
2348         const ZZ_pEXTransMultiplier& B, const ZZ_pEXModulus& F)
2349{
2350   ZZ_pEX xx;
2351   TransMulMod(xx, to_ZZ_pEX(a), B, F);
2352   x = xx.rep;
2353}
2354
2355static
2356void ProjectPowers(vec_ZZ_pE& x, const ZZ_pEX& a, long k, 
2357                   const ZZ_pEXArgument& H, const ZZ_pEXModulus& F)
2358{
2359   if (k < 0 || NTL_OVERFLOW(k, 1, 0) || deg(a) >= F.n) 
2360      Error("ProjectPowers: bad args");
2361
2362   long m = H.H.length()-1;
2363   long l = (k+m-1)/m - 1;
2364
2365   ZZ_pEXTransMultiplier M;
2366   build(M, H.H[m], F);
2367
2368   ZZ_pEX s;
2369   s = a;
2370
2371   x.SetLength(k);
2372
2373   long i;
2374
2375   for (i = 0; i <= l; i++) {
2376      long m1 = min(m, k-i*m);
2377      for (long j = 0; j < m1; j++)
2378         InnerProduct(x[i*m+j], H.H[j].rep, s.rep);
2379      if (i < l)
2380         TransMulMod(s, s, M, F);
2381   }
2382}
2383
2384static
2385void ProjectPowers(vec_ZZ_pE& x, const ZZ_pEX& a, long k, const ZZ_pEX& h, 
2386                   const ZZ_pEXModulus& F)
2387{
2388   if (k < 0 || deg(a) >= F.n || deg(h) >= F.n)
2389      Error("ProjectPowers: bad args");
2390
2391   if (k == 0) {
2392      x.SetLength(0);;
2393      return;
2394   }
2395
2396   long m = SqrRoot(k);
2397
2398   ZZ_pEXArgument H;
2399   build(H, h, F, m);
2400
2401   ProjectPowers(x, a, k, H, F);
2402}
2403
2404void ProjectPowers(vec_ZZ_pE& x, const vec_ZZ_pE& a, long k,
2405                   const ZZ_pEXArgument& H, const ZZ_pEXModulus& F)
2406{
2407   ProjectPowers(x, to_ZZ_pEX(a), k, H, F);
2408}
2409
2410void ProjectPowers(vec_ZZ_pE& x, const vec_ZZ_pE& a, long k,
2411                   const ZZ_pEX& h, const ZZ_pEXModulus& F)
2412{
2413   ProjectPowers(x, to_ZZ_pEX(a), k, h, F);
2414}
2415
2416
2417
2418
2419void BerlekampMassey(ZZ_pEX& h, const vec_ZZ_pE& a, long m)
2420{
2421   ZZ_pEX Lambda, Sigma, Temp;
2422   long L;
2423   ZZ_pE Delta, Delta1, t1;
2424   long shamt;
2425
2426   // cerr << "*** " << m << "\n";
2427
2428   Lambda.SetMaxLength(m+1);
2429   Sigma.SetMaxLength(m+1);
2430   Temp.SetMaxLength(m+1);
2431
2432   L = 0;
2433   set(Lambda);
2434   clear(Sigma);
2435   set(Delta);
2436   shamt = 0;
2437
2438   long i, r, dl;
2439
2440   for (r = 1; r <= 2*m; r++) {
2441      // cerr << r << "--";
2442      clear(Delta1);
2443      dl = deg(Lambda);
2444      for (i = 0; i <= dl; i++) {
2445         mul(t1, Lambda.rep[i], a[r-i-1]);
2446         add(Delta1, Delta1, t1);
2447      }
2448
2449      if (IsZero(Delta1)) {
2450         shamt++;
2451         // cerr << "case 1: " << deg(Lambda) << " " << deg(Sigma) << " " << shamt << "\n";
2452      }
2453      else if (2*L < r) {
2454         div(t1, Delta1, Delta);
2455         mul(Temp, Sigma, t1);
2456         Sigma = Lambda;
2457         ShiftSub(Lambda, Temp, shamt+1);
2458         shamt = 0;
2459         L = r-L;
2460         Delta = Delta1;
2461         // cerr << "case 2: " << deg(Lambda) << " " << deg(Sigma) << " " << shamt << "\n";
2462      }
2463      else {
2464         shamt++;
2465         div(t1, Delta1, Delta);
2466         mul(Temp, Sigma, t1);
2467         ShiftSub(Lambda, Temp, shamt);
2468         // cerr << "case 3: " << deg(Lambda) << " " << deg(Sigma) << " " << shamt << "\n";
2469      }
2470   }
2471
2472   // cerr << "finished: " << L << " " << deg(Lambda) << "\n";
2473
2474   dl = deg(Lambda);
2475   h.rep.SetLength(L + 1);
2476
2477   for (i = 0; i < L - dl; i++)
2478      clear(h.rep[i]);
2479
2480   for (i = L - dl; i <= L; i++)
2481      h.rep[i] = Lambda.rep[L - i];
2482}
2483
2484
2485
2486
2487void MinPolySeq(ZZ_pEX& h, const vec_ZZ_pE& a, long m)
2488{
2489   if (m < 0 || NTL_OVERFLOW(m, 1, 0)) Error("MinPoly: bad args");
2490   if (a.length() < 2*m) Error("MinPoly: sequence too short");
2491
2492   BerlekampMassey(h, a, m);
2493}
2494
2495
2496void DoMinPolyMod(ZZ_pEX& h, const ZZ_pEX& g, const ZZ_pEXModulus& F, long m, 
2497               const ZZ_pEX& R)
2498{
2499   vec_ZZ_pE x;
2500
2501   ProjectPowers(x, R, 2*m, g, F);
2502   MinPolySeq(h, x, m);
2503}
2504
2505void ProbMinPolyMod(ZZ_pEX& h, const ZZ_pEX& g, const ZZ_pEXModulus& F, long m)
2506{
2507   long n = F.n;
2508   if (m < 1 || m > n) Error("ProbMinPoly: bad args");
2509
2510   ZZ_pEX R;
2511   random(R, n);
2512
2513   DoMinPolyMod(h, g, F, m, R);
2514}
2515
2516void ProbMinPolyMod(ZZ_pEX& h, const ZZ_pEX& g, const ZZ_pEXModulus& F)
2517{
2518   ProbMinPolyMod(h, g, F, F.n);
2519}
2520
2521void MinPolyMod(ZZ_pEX& hh, const ZZ_pEX& g, const ZZ_pEXModulus& F, long m)
2522{
2523   ZZ_pEX h, h1;
2524   long n = F.n;
2525   if (m < 1 || m > n) Error("MinPoly: bad args");
2526
2527   /* probabilistically compute min-poly */
2528
2529   ProbMinPolyMod(h, g, F, m);
2530   if (deg(h) == m) { hh = h; return; }
2531   CompMod(h1, h, g, F);
2532   if (IsZero(h1)) { hh = h; return; }
2533
2534   /* not completely successful...must iterate */
2535
2536   ZZ_pEX h2, h3;
2537   ZZ_pEX R;
2538   ZZ_pEXTransMultiplier H1;
2539   
2540
2541   for (;;) {
2542      random(R, n);
2543      build(H1, h1, F);
2544      TransMulMod(R, R, H1, F);
2545      DoMinPolyMod(h2, g, F, m-deg(h), R);
2546
2547      mul(h, h, h2);
2548      if (deg(h) == m) { hh = h; return; }
2549      CompMod(h3, h2, g, F);
2550      MulMod(h1, h3, h1, F);
2551      if (IsZero(h1)) { hh = h; return; }
2552   }
2553}
2554
2555void IrredPolyMod(ZZ_pEX& h, const ZZ_pEX& g, const ZZ_pEXModulus& F, long m)
2556{
2557   if (m < 1 || m > F.n) Error("IrredPoly: bad args");
2558
2559   ZZ_pEX R;
2560   set(R);
2561
2562   DoMinPolyMod(h, g, F, m, R);
2563}
2564
2565
2566
2567void IrredPolyMod(ZZ_pEX& h, const ZZ_pEX& g, const ZZ_pEXModulus& F)
2568{
2569   IrredPolyMod(h, g, F, F.n);
2570}
2571
2572
2573
2574void MinPolyMod(ZZ_pEX& hh, const ZZ_pEX& g, const ZZ_pEXModulus& F)
2575{
2576   MinPolyMod(hh, g, F, F.n);
2577}
2578
2579void diff(ZZ_pEX& x, const ZZ_pEX& a)
2580{
2581   long n = deg(a);
2582   long i;
2583
2584   if (n <= 0) {
2585      clear(x);
2586      return;
2587   }
2588
2589   if (&x != &a)
2590      x.rep.SetLength(n);
2591
2592   for (i = 0; i <= n-1; i++) {
2593      mul(x.rep[i], a.rep[i+1], i+1);
2594   }
2595
2596   if (&x == &a)
2597      x.rep.SetLength(n);
2598
2599   x.normalize();
2600}
2601
2602
2603
2604void MakeMonic(ZZ_pEX& x)
2605{
2606   if (IsZero(x))
2607      return;
2608
2609   if (IsOne(LeadCoeff(x)))
2610      return;
2611
2612   ZZ_pE t;
2613
2614   inv(t, LeadCoeff(x));
2615   mul(x, x, t);
2616}
2617
2618
2619long divide(ZZ_pEX& q, const ZZ_pEX& a, const ZZ_pEX& b)
2620{
2621   if (IsZero(b)) {
2622      if (IsZero(a)) {
2623         clear(q);
2624         return 1;
2625      }
2626      else
2627         return 0;
2628   }
2629
2630   ZZ_pEX lq, r;
2631   DivRem(lq, r, a, b);
2632   if (!IsZero(r)) return 0; 
2633   q = lq;
2634   return 1;
2635}
2636
2637long divide(const ZZ_pEX& a, const ZZ_pEX& b)
2638{
2639   if (IsZero(b)) return IsZero(a);
2640   ZZ_pEX lq, r;
2641   DivRem(lq, r, a, b);
2642   if (!IsZero(r)) return 0; 
2643   return 1;
2644}
2645
2646
2647
2648static
2649long OptWinSize(long n)
2650// finds k that minimizes n/(k+1) + 2^{k-1}
2651
2652{
2653   long k;
2654   double v, v_new;
2655
2656
2657   v = n/2.0 + 1.0;
2658   k = 1;
2659
2660   for (;;) {
2661      v_new = n/(double(k+2)) + double(1L << k);
2662      if (v_new >= v) break;
2663      v = v_new;
2664      k++;
2665   }
2666
2667   return k;
2668}
2669     
2670
2671
2672void PowerMod(ZZ_pEX& h, const ZZ_pEX& g, const ZZ& e, const ZZ_pEXModulus& F)
2673// h = g^e mod f using "sliding window" algorithm
2674{
2675   if (deg(g) >= F.n) Error("PowerMod: bad args");
2676
2677   if (e == 0) {
2678      set(h);
2679      return;
2680   }
2681
2682   if (e == 1) {
2683      h = g;
2684      return;
2685   }
2686
2687   if (e == -1) {
2688      InvMod(h, g, F);
2689      return;
2690   }
2691
2692   if (e == 2) {
2693      SqrMod(h, g, F);
2694      return;
2695   }
2696
2697   if (e == -2) {
2698      SqrMod(h, g, F);
2699      InvMod(h, h, F);
2700      return;
2701   }
2702
2703
2704   long n = NumBits(e);
2705
2706   ZZ_pEX res;
2707   res.SetMaxLength(F.n);
2708   set(res);
2709
2710   long i;
2711
2712   if (n < 16) {
2713      // plain square-and-multiply algorithm
2714
2715      for (i = n - 1; i >= 0; i--) {
2716         SqrMod(res, res, F);
2717         if (bit(e, i))
2718            MulMod(res, res, g, F);
2719      }
2720
2721      if (e < 0) InvMod(res, res, F);
2722
2723      h = res;
2724      return;
2725   }
2726
2727   long k = OptWinSize(n);
2728   k = min(k, 3);
2729
2730   vec_ZZ_pEX v;
2731
2732   v.SetLength(1L << (k-1));
2733
2734   v[0] = g;
2735 
2736   if (k > 1) {
2737      ZZ_pEX t;
2738      SqrMod(t, g, F);
2739
2740      for (i = 1; i < (1L << (k-1)); i++)
2741         MulMod(v[i], v[i-1], t, F);
2742   }
2743
2744
2745   long val;
2746   long cnt;
2747   long m;
2748
2749   val = 0;
2750   for (i = n-1; i >= 0; i--) {
2751      val = (val << 1) | bit(e, i); 
2752      if (val == 0)
2753         SqrMod(res, res, F);
2754      else if (val >= (1L << (k-1)) || i == 0) {
2755         cnt = 0;
2756         while ((val & 1) == 0) {
2757            val = val >> 1;
2758            cnt++;
2759         }
2760
2761         m = val;
2762         while (m > 0) {
2763            SqrMod(res, res, F);
2764            m = m >> 1;
2765         }
2766
2767         MulMod(res, res, v[val >> 1], F);
2768
2769         while (cnt > 0) {
2770            SqrMod(res, res, F);
2771            cnt--;
2772         }
2773
2774         val = 0;
2775      }
2776   }
2777
2778   if (e < 0) InvMod(res, res, F);
2779
2780   h = res;
2781}
2782
2783void InvMod(ZZ_pEX& x, const ZZ_pEX& a, const ZZ_pEX& f)
2784{
2785   if (deg(a) >= deg(f) || deg(f) == 0) Error("InvMod: bad args");
2786
2787   ZZ_pEX d, t;
2788
2789   XGCD(d, x, t, a, f);
2790   if (!IsOne(d))
2791      Error("ZZ_pEX InvMod: can't compute multiplicative inverse");
2792}
2793
2794long InvModStatus(ZZ_pEX& x, const ZZ_pEX& a, const ZZ_pEX& f)
2795{
2796   if (deg(a) >= deg(f) || deg(f) == 0) Error("InvModStatus: bad args");
2797   ZZ_pEX d, t;
2798
2799   XGCD(d, x, t, a, f);
2800   if (!IsOne(d)) {
2801      x = d;
2802      return 1;
2803   }
2804   else
2805      return 0;
2806}
2807
2808
2809void MulMod(ZZ_pEX& x, const ZZ_pEX& a, const ZZ_pEX& b, const ZZ_pEX& f)
2810{
2811   if (deg(a) >= deg(f) || deg(b) >= deg(f) || deg(f) == 0)
2812      Error("MulMod: bad args");
2813
2814   ZZ_pEX t;
2815
2816   mul(t, a, b);
2817   rem(x, t, f);
2818}
2819
2820void SqrMod(ZZ_pEX& x, const ZZ_pEX& a, const ZZ_pEX& f)
2821{
2822   if (deg(a) >= deg(f) || deg(f) == 0) Error("SqrMod: bad args");
2823
2824   ZZ_pEX t;
2825
2826   sqr(t, a);
2827   rem(x, t, f);
2828}
2829
2830
2831void PowerXMod(ZZ_pEX& hh, const ZZ& e, const ZZ_pEXModulus& F)
2832{
2833   if (F.n < 0) Error("PowerXMod: uninitialized modulus");
2834
2835   if (IsZero(e)) {
2836      set(hh);
2837      return;
2838   }
2839
2840   long n = NumBits(e);
2841   long i;
2842
2843   ZZ_pEX h;
2844
2845   h.SetMaxLength(F.n);
2846   set(h);
2847
2848   for (i = n - 1; i >= 0; i--) {
2849      SqrMod(h, h, F);
2850      if (bit(e, i))
2851         MulByXMod(h, h, F.f);
2852   }
2853
2854   if (e < 0) InvMod(h, h, F);
2855
2856   hh = h;
2857}
2858
2859
2860void reverse(ZZ_pEX& x, const ZZ_pEX& a, long hi)
2861{
2862   if (hi < 0) { clear(x); return; }
2863   if (NTL_OVERFLOW(hi, 1, 0))
2864      Error("overflow in reverse");
2865
2866   if (&x == &a) {
2867      ZZ_pEX tmp;
2868      CopyReverse(tmp, a, hi);
2869      x = tmp;
2870   }
2871   else
2872      CopyReverse(x, a, hi);
2873}
2874
2875
2876void power(ZZ_pEX& x, const ZZ_pEX& a, long e)
2877{
2878   if (e < 0) {
2879      Error("power: negative exponent");
2880   }
2881
2882   if (e == 0) {
2883      x = 1;
2884      return;
2885   }
2886
2887   if (a == 0 || a == 1) {
2888      x = a;
2889      return;
2890   }
2891
2892   long da = deg(a);
2893
2894   if (da == 0) {
2895      x = power(ConstTerm(a), e);
2896      return;
2897   }
2898
2899   if (da > (NTL_MAX_LONG-1)/e)
2900      Error("overflow in power");
2901
2902   ZZ_pEX res;
2903   res.SetMaxLength(da*e + 1);
2904   res = 1;
2905   
2906   long k = NumBits(e);
2907   long i;
2908
2909   for (i = k - 1; i >= 0; i--) {
2910      sqr(res, res);
2911      if (bit(e, i))
2912         mul(res, res, a);
2913   }
2914
2915   x = res;
2916}
2917
2918
2919
2920static
2921void FastTraceVec(vec_ZZ_pE& S, const ZZ_pEXModulus& f)
2922{
2923   long n = deg(f);
2924
2925   ZZ_pEX x = reverse(-LeftShift(reverse(diff(reverse(f)), n-1), n-1)/f, n-1);
2926
2927   S.SetLength(n);
2928   S[0] = n;
2929
2930   long i;
2931   for (i = 1; i < n; i++)
2932      S[i] = coeff(x, i);
2933}
2934
2935
2936void PlainTraceVec(vec_ZZ_pE& S, const ZZ_pEX& ff)
2937{
2938   if (deg(ff) <= 0)
2939      Error("TraceVec: bad args");
2940
2941   ZZ_pEX f;
2942   f = ff;
2943
2944   MakeMonic(f);
2945
2946   long n = deg(f);
2947
2948   S.SetLength(n);
2949
2950   if (n == 0)
2951      return;
2952
2953   long k, i;
2954   ZZ_pX acc, t;
2955   ZZ_pE t1;
2956
2957   S[0] = n;
2958
2959   for (k = 1; k < n; k++) {
2960      mul(acc, rep(f.rep[n-k]), k);
2961
2962      for (i = 1; i < k; i++) {
2963         mul(t, rep(f.rep[n-i]), rep(S[k-i]));
2964         add(acc, acc, t);
2965      }
2966
2967      conv(t1, acc);
2968      negate(S[k], t1);
2969   }
2970}
2971
2972void TraceVec(vec_ZZ_pE& S, const ZZ_pEX& f)
2973{
2974   if (deg(f) < ZZ_pE::DivCross())
2975      PlainTraceVec(S, f);
2976   else
2977      FastTraceVec(S, f);
2978}
2979
2980static
2981void ComputeTraceVec(const ZZ_pEXModulus& F)
2982{
2983   vec_ZZ_pE& S = *((vec_ZZ_pE *) &F.tracevec);
2984
2985   if (S.length() > 0)
2986      return;
2987
2988   if (F.method == ZZ_pEX_MOD_PLAIN) {
2989      PlainTraceVec(S, F.f);
2990   }
2991   else {
2992      FastTraceVec(S, F);
2993   }
2994}
2995
2996void TraceMod(ZZ_pE& x, const ZZ_pEX& a, const ZZ_pEXModulus& F)
2997{
2998   long n = F.n;
2999
3000   if (deg(a) >= n)
3001      Error("trace: bad args");
3002
3003   if (F.tracevec.length() == 0) 
3004      ComputeTraceVec(F);
3005
3006   InnerProduct(x, a.rep, F.tracevec);
3007}
3008
3009void TraceMod(ZZ_pE& x, const ZZ_pEX& a, const ZZ_pEX& f)
3010{
3011   if (deg(a) >= deg(f) || deg(f) <= 0)
3012      Error("trace: bad args");
3013
3014   project(x, TraceVec(f), a);
3015}
3016
3017
3018void PlainResultant(ZZ_pE& rres, const ZZ_pEX& a, const ZZ_pEX& b)
3019{
3020   ZZ_pE res;
3021 
3022   if (IsZero(a) || IsZero(b))
3023      clear(res);
3024   else if (deg(a) == 0 && deg(b) == 0) 
3025      set(res);
3026   else {
3027      long d0, d1, d2;
3028      ZZ_pE lc;
3029      set(res);
3030
3031      long n = max(deg(a),deg(b)) + 1;
3032      ZZ_pEX u(INIT_SIZE, n), v(INIT_SIZE, n);
3033      vec_ZZ_pX tmp;
3034      SetSize(tmp, n, 2*ZZ_pE::degree());
3035
3036      u = a;
3037      v = b;
3038
3039      for (;;) {
3040         d0 = deg(u);
3041         d1 = deg(v);
3042         lc = LeadCoeff(v);
3043
3044         PlainRem(u, u, v, tmp);
3045         swap(u, v);
3046
3047         d2 = deg(v);
3048         if (d2 >= 0) {
3049            power(lc, lc, d0-d2);
3050            mul(res, res, lc);
3051            if (d0 & d1 & 1) negate(res, res);
3052         }
3053         else {
3054            if (d1 == 0) {
3055               power(lc, lc, d0);
3056               mul(res, res, lc);
3057            }
3058            else
3059               clear(res);
3060       
3061            break;
3062         }
3063      }
3064
3065      rres = res;
3066   }
3067}
3068
3069void resultant(ZZ_pE& rres, const ZZ_pEX& a, const ZZ_pEX& b)
3070{
3071   PlainResultant(rres, a, b); 
3072}
3073
3074
3075void NormMod(ZZ_pE& x, const ZZ_pEX& a, const ZZ_pEX& f)
3076{
3077   if (deg(f) <= 0 || deg(a) >= deg(f)) 
3078      Error("norm: bad args");
3079
3080   if (IsZero(a)) {
3081      clear(x);
3082      return;
3083   }
3084
3085   ZZ_pE t;
3086   resultant(t, f, a);
3087   if (!IsOne(LeadCoeff(f))) {
3088      ZZ_pE t1;
3089      power(t1, LeadCoeff(f), deg(a));
3090      inv(t1, t1);
3091      mul(t, t, t1);
3092   }
3093
3094   x = t;
3095}
3096
3097
3098
3099// tower stuff...
3100
3101
3102
3103void InnerProduct(ZZ_pEX& x, const vec_ZZ_p& v, long low, long high,
3104                   const vec_ZZ_pEX& H, long n, vec_ZZ_pE& t)
3105{
3106   ZZ_pE s;
3107   long i, j;
3108
3109   for (j = 0; j < n; j++)
3110      clear(t[j]);
3111
3112   high = min(high, v.length()-1);
3113   for (i = low; i <= high; i++) {
3114      const vec_ZZ_pE& h = H[i-low].rep;
3115      long m = h.length();
3116      const ZZ_p& w = v[i];
3117
3118      for (j = 0; j < m; j++) {
3119         mul(s, h[j], w);
3120         add(t[j], t[j], s);
3121      }
3122   }
3123
3124   x.rep.SetLength(n);
3125   for (j = 0; j < n; j++)
3126      x.rep[j] = t[j];
3127
3128   x.normalize();
3129}
3130
3131
3132
3133void CompTower(ZZ_pEX& x, const ZZ_pX& g, const ZZ_pEXArgument& A,
3134             const ZZ_pEXModulus& F)
3135{
3136   if (deg(g) <= 0) {
3137      conv(x, g);
3138      return;
3139   }
3140
3141
3142   ZZ_pEX s, t;
3143   vec_ZZ_pE scratch;
3144   scratch.SetLength(deg(F));
3145
3146   long m = A.H.length() - 1;
3147   long l = ((g.rep.length()+m-1)/m) - 1;
3148
3149   const ZZ_pEX& M = A.H[m];
3150
3151   InnerProduct(t, g.rep, l*m, l*m + m - 1, A.H, F.n, scratch);
3152   for (long i = l-1; i >= 0; i--) {
3153      InnerProduct(s, g.rep, i*m, i*m + m - 1, A.H, F.n, scratch);
3154      MulMod(t, t, M, F);
3155      add(t, t, s);
3156   }
3157   x = t;
3158}
3159
3160
3161void CompTower(ZZ_pEX& x, const ZZ_pX& g, const ZZ_pEX& h, 
3162             const ZZ_pEXModulus& F)
3163   // x = g(h) mod f
3164{
3165   long m = SqrRoot(g.rep.length());
3166
3167   if (m == 0) {
3168      clear(x);
3169      return;
3170   }
3171
3172
3173   ZZ_pEXArgument A;
3174
3175   build(A, h, F, m);
3176
3177   CompTower(x, g, A, F);
3178}
3179
3180void PrepareProjection(vec_vec_ZZ_p& tt, const vec_ZZ_pE& s,
3181                       const vec_ZZ_p& proj)
3182{
3183   long l = s.length();
3184   tt.SetLength(l);
3185
3186   ZZ_pXMultiplier M;
3187   long i;
3188
3189   for (i = 0; i < l; i++) {
3190      build(M, rep(s[i]), ZZ_pE::modulus());
3191      UpdateMap(tt[i], proj, M, ZZ_pE::modulus());
3192   }
3193}
3194
3195void ProjectedInnerProduct(ZZ_p& x, const vec_ZZ_pE& a, 
3196                           const vec_vec_ZZ_p& b)
3197{
3198   long n = min(a.length(), b.length());
3199
3200   ZZ_p t, res;
3201
3202   res = 0;
3203
3204   long i;
3205   for (i = 0; i < n; i++) {
3206      project(t, b[i], rep(a[i]));
3207      res += t;
3208   }
3209
3210   x = res;
3211}
3212
3213
3214void PrecomputeProj(vec_ZZ_p& proj, const ZZ_pX& f)
3215{
3216   long n = deg(f);
3217
3218   if (n <= 0) Error("PrecomputeProj: bad args");
3219
3220   if (ConstTerm(f) != 0) {
3221      proj.SetLength(1);
3222      proj[0] = 1;
3223   }
3224   else {
3225      proj.SetLength(n);
3226      clear(proj);
3227      proj[n-1] = 1;
3228   }
3229}
3230   
3231
3232
3233void ProjectPowersTower(vec_ZZ_p& x, const vec_ZZ_pE& a, long k,
3234                   const ZZ_pEXArgument& H, const ZZ_pEXModulus& F,
3235                   const vec_ZZ_p& proj)
3236
3237{
3238   long n = F.n;
3239
3240   if (a.length() > n || k < 0 || NTL_OVERFLOW(k, 1, 0)) 
3241      Error("ProjectPowers: bad args");
3242
3243   long m = H.H.length()-1;
3244   long l = (k+m-1)/m - 1;
3245
3246   ZZ_pEXTransMultiplier M;
3247   build(M, H.H[m], F);
3248
3249   vec_ZZ_pE s(INIT_SIZE, n);
3250   s = a;
3251
3252   x.SetLength(k);
3253
3254   vec_vec_ZZ_p tt;
3255
3256   for (long i = 0; i <= l; i++) {
3257      long m1 = min(m, k-i*m);
3258      ZZ_p* w = &x[i*m];
3259
3260      PrepareProjection(tt, s, proj);
3261
3262      for (long j = 0; j < m1; j++)
3263         ProjectedInnerProduct(w[j], H.H[j].rep, tt);
3264      if (i < l)
3265         UpdateMap(s, s, M, F);
3266   }
3267}
3268
3269
3270
3271
3272void ProjectPowersTower(vec_ZZ_p& x, const vec_ZZ_pE& a, long k,
3273                   const ZZ_pEX& h, const ZZ_pEXModulus& F,
3274                   const vec_ZZ_p& proj)
3275
3276{
3277   if (a.length() > F.n || k < 0) Error("ProjectPowers: bad args");
3278
3279   if (k == 0) {
3280      x.SetLength(0);
3281      return;
3282   }
3283
3284   long m = SqrRoot(k);
3285
3286   ZZ_pEXArgument H;
3287
3288   build(H, h, F, m);
3289   ProjectPowersTower(x, a, k, H, F, proj);
3290}
3291
3292
3293void DoMinPolyTower(ZZ_pX& h, const ZZ_pEX& g, const ZZ_pEXModulus& F, long m,
3294               const vec_ZZ_pE& R, const vec_ZZ_p& proj)
3295{
3296   vec_ZZ_p x;
3297
3298   ProjectPowersTower(x, R, 2*m, g, F, proj);
3299   
3300   MinPolySeq(h, x, m);
3301}
3302
3303
3304void ProbMinPolyTower(ZZ_pX& h, const ZZ_pEX& g, const ZZ_pEXModulus& F, 
3305                      long m)
3306{
3307   long n = F.n;
3308   if (m < 1 || m > n*ZZ_pE::degree())
3309      Error("MinPoly: bad args");
3310
3311   vec_ZZ_pE R;
3312   R.SetLength(n);
3313   long i;
3314   for (i = 0; i < n; i++)
3315      random(R[i]);
3316
3317   vec_ZZ_p proj;
3318   PrecomputeProj(proj, ZZ_pE::modulus());
3319
3320   DoMinPolyTower(h, g, F, m, R, proj);
3321}
3322
3323
3324void ProbMinPolyTower(ZZ_pX& h, const ZZ_pEX& g, const ZZ_pEXModulus& F, 
3325                      long m, const vec_ZZ_p& proj)
3326{
3327   long n = F.n;
3328   if (m < 1 || m > n*ZZ_pE::degree())
3329      Error("MinPoly: bad args");
3330
3331   vec_ZZ_pE R;
3332   R.SetLength(n);
3333   long i;
3334   for (i = 0; i < n; i++)
3335      random(R[i]);
3336
3337   DoMinPolyTower(h, g, F, m, R, proj);
3338}
3339
3340void MinPolyTower(ZZ_pX& hh, const ZZ_pEX& g, const ZZ_pEXModulus& F, long m)
3341{
3342   ZZ_pX h;
3343   ZZ_pEX h1;
3344   long n = F.n;
3345
3346   if (m < 1 || m > n*ZZ_pE::degree())
3347      Error("MinPoly: bad args");
3348
3349   vec_ZZ_p proj;
3350   PrecomputeProj(proj, ZZ_pE::modulus());
3351
3352   /* probabilistically compute min-poly */
3353
3354   ProbMinPolyTower(h, g, F, m, proj);
3355   if (deg(h) == m) { hh = h; return; }
3356   CompTower(h1, h, g, F);
3357   if (IsZero(h1)) { hh = h; return; }
3358
3359   /* not completely successful...must iterate */
3360
3361   long i;
3362
3363   ZZ_pX h2;
3364   ZZ_pEX h3;
3365   vec_ZZ_pE R;
3366   ZZ_pEXTransMultiplier H1;
3367   
3368
3369   for (;;) {
3370      R.SetLength(n);
3371      for (i = 0; i < n; i++) random(R[i]);
3372      build(H1, h1, F);
3373      UpdateMap(R, R, H1, F);
3374      DoMinPolyTower(h2, g, F, m-deg(h), R, proj);
3375
3376      mul(h, h, h2);
3377      if (deg(h) == m) { hh = h; return; }
3378      CompTower(h3, h2, g, F);
3379      MulMod(h1, h3, h1, F);
3380      if (IsZero(h1)) { hh = h; return; }
3381   }
3382}
3383
3384void IrredPolyTower(ZZ_pX& h, const ZZ_pEX& g, const ZZ_pEXModulus& F, long m)
3385{
3386   if (m < 1 || m > deg(F)*ZZ_pE::degree())
3387      Error("IrredPoly: bad args");
3388
3389   vec_ZZ_pE R;
3390   R.SetLength(1);
3391   R[0] = 1;
3392
3393   vec_ZZ_p proj;
3394   proj.SetLength(1);
3395   proj[0] = 1;
3396
3397   DoMinPolyTower(h, g, F, m, R, proj);
3398}
3399
3400NTL_END_IMPL
Note: See TracBrowser for help on using the repository browser.