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

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