source: git/ntl/src/ZZ_pX.c @ de6a29

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