source: git/ntl/src/lzz_pX.c @ 33a041

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