source: git/ntl/src/ZZ_pEX.c @ 26e030

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