source: git/ntl/src/lzz_pEX.c @ eb5966

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