source: git/ntl/src/lzz_pEX.c @ 45fefa

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