source: git/ntl/src/lzz_pEX.c @ 2cfffe

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