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

spielwiese
Last change on this file since de6a29 was de6a29, checked in by Hans Schönemann <hannes@…>, 18 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: 34.4 KB
Line 
1
2
3#include <NTL/lzz_pX.h>
4
5#include <NTL/new.h>
6
7NTL_START_IMPL
8
9
10
11long divide(zz_pX& q, const zz_pX& a, const zz_pX& b)
12{
13   if (IsZero(b)) {
14      if (IsZero(a)) {
15         clear(q);
16         return 1;
17      }
18      else
19         return 0;
20   }
21
22   zz_pX lq, r;
23   DivRem(lq, r, a, b);
24   if (!IsZero(r)) return 0;
25   q = lq;
26   return 1;
27}
28
29long divide(const zz_pX& a, const zz_pX& b)
30{
31   if (IsZero(b)) return IsZero(a);
32   zz_pX lq, r;
33   DivRem(lq, r, a, b);
34   if (!IsZero(r)) return 0;
35   return 1;
36}
37
38
39
40void zz_pXMatrix::operator=(const zz_pXMatrix& M)
41{
42   elts[0][0] = M.elts[0][0];
43   elts[0][1] = M.elts[0][1];
44   elts[1][0] = M.elts[1][0];
45   elts[1][1] = M.elts[1][1];
46}
47
48
49void RightShift(zz_pX& x, const zz_pX& a, long n)
50{
51   if (IsZero(a)) {
52      clear(x);
53      return;
54   }
55
56   if (n < 0) {
57      if (n < -NTL_MAX_LONG) Error("overflow in RightShift");
58      LeftShift(x, a, -n);
59      return;
60   }
61
62   long da = deg(a);
63   long i;
64 
65   if (da < n) {
66      clear(x);
67      return;
68   }
69
70   if (&x != &a)
71      x.rep.SetLength(da-n+1);
72
73   for (i = 0; i <= da-n; i++)
74      x.rep[i] = a.rep[i+n];
75
76   if (&x == &a)
77      x.rep.SetLength(da-n+1);
78
79   x.normalize();
80}
81
82void LeftShift(zz_pX& x, const zz_pX& a, long n)
83{
84   if (IsZero(a)) {
85      clear(x);
86      return;
87   }
88
89   if (n < 0) {
90      if (n < -NTL_MAX_LONG) 
91         clear(x);
92      else
93         RightShift(x, a, -n);
94      return;
95   }
96
97   if (NTL_OVERFLOW(n, 1, 0))
98      Error("overflow in LeftShift");
99
100   long m = a.rep.length();
101
102   x.rep.SetLength(m+n);
103
104   long i;
105   for (i = m-1; i >= 0; i--)
106      x.rep[i+n] = a.rep[i];
107
108   for (i = 0; i < n; i++)
109      clear(x.rep[i]);
110}
111
112
113void ShiftAdd(zz_pX& U, const zz_pX& V, long n)
114// assumes input does not alias output
115{
116   if (IsZero(V))
117      return;
118
119   long du = deg(U);
120   long dv = deg(V);
121
122   long d = max(du, n+dv);
123
124   U.rep.SetLength(d+1);
125   long i;
126
127   for (i = du+1; i <= d; i++)
128      clear(U.rep[i]);
129
130   for (i = 0; i <= dv; i++)
131      add(U.rep[i+n], U.rep[i+n], V.rep[i]);
132
133   U.normalize();
134}
135
136void ShiftSub(zz_pX& U, const zz_pX& V, long n)
137// assumes input does not alias output
138{
139   if (IsZero(V))
140      return;
141
142   long du = deg(U);
143   long dv = deg(V);
144
145   long d = max(du, n+dv);
146
147   U.rep.SetLength(d+1);
148   long i;
149
150   for (i = du+1; i <= d; i++)
151      clear(U.rep[i]);
152
153   for (i = 0; i <= dv; i++)
154      sub(U.rep[i+n], U.rep[i+n], V.rep[i]);
155
156   U.normalize();
157}
158
159void mul(zz_pX& U, zz_pX& V, const zz_pXMatrix& M)
160// (U, V)^T = M*(U, V)^T
161{
162   long d = deg(U) - deg(M(1,1));
163   long k = NextPowerOfTwo(d - 1);
164
165   // When the GCD algorithm is run on polynomials of degree n, n-1,
166   // where n is a power of two, then d-1 is likely to be a power of two.
167   // It would be more natural to set k = NextPowerOfTwo(d+1), but this
168   // would be much less efficient in this case.
169
170   long n = (1L << k);
171   long xx;
172   zz_p a0, a1, b0, b1, c0, d0, u0, u1, v0, v1, nu0, nu1, nv0;
173   zz_p t1, t2;
174
175   if (n == d-1)
176      xx = 1;
177   else if (n == d)
178      xx = 2;
179   else 
180      xx = 3;
181
182   switch (xx) {
183   case 1:
184      GetCoeff(a0, M(0,0), 0);
185      GetCoeff(a1, M(0,0), 1);
186      GetCoeff(b0, M(0,1), 0);
187      GetCoeff(b1, M(0,1), 1);
188      GetCoeff(c0, M(1,0), 0);
189      GetCoeff(d0, M(1,1), 0);
190
191      GetCoeff(u0, U, 0);
192      GetCoeff(u1, U, 1);
193      GetCoeff(v0, V, 0);
194      GetCoeff(v1, V, 1);
195
196      mul(t1, (a0), (u0));
197      mul(t2, (b0), (v0));
198      add(t1, t1, t2); 
199      nu0 = t1;
200
201      mul(t1, (a1), (u0));
202      mul(t2, (a0), (u1));
203      add(t1, t1, t2);
204      mul(t2, (b1), (v0));
205      add(t1, t1, t2);
206      mul(t2, (b0), (v1));
207      add(t1, t1, t2);
208      nu1 = t1;
209
210      mul(t1, (c0), (u0));
211      mul(t2, (d0), (v0));
212      add (t1, t1, t2);
213      nv0 = t1;
214   
215      break;
216
217   case 2:
218      GetCoeff(a0, M(0,0), 0);
219      GetCoeff(b0, M(0,1), 0);
220
221      GetCoeff(u0, U, 0);
222      GetCoeff(v0, V, 0);
223
224      mul(t1, (a0), (u0));
225      mul(t2, (b0), (v0));
226      add(t1, t1, t2); 
227      nu0 = t1;
228
229      break;
230
231   case 3:
232      break;
233
234   }
235
236   fftRep RU(INIT_SIZE, k), RV(INIT_SIZE, k), R1(INIT_SIZE, k), 
237          R2(INIT_SIZE, k);
238
239   TofftRep(RU, U, k); 
240   TofftRep(RV, V, k); 
241
242   TofftRep(R1, M(0,0), k);
243   mul(R1, R1, RU);
244   TofftRep(R2, M(0,1), k);
245   mul(R2, R2, RV);
246   add(R1, R1, R2);
247   FromfftRep(U, R1, 0, d);
248
249   TofftRep(R1, M(1,0), k);
250   mul(R1, R1, RU);
251   TofftRep(R2, M(1,1), k);
252   mul(R2, R2, RV);
253   add(R1, R1, R2);
254   FromfftRep(V, R1, 0, d-1);
255
256   // now fix-up results
257
258   switch (xx) {
259   case 1:
260      GetCoeff(u0, U, 0);
261      sub(u0, u0, nu0);
262      SetCoeff(U, d-1, u0);
263      SetCoeff(U, 0, nu0);
264
265      GetCoeff(u1, U, 1);
266      sub(u1, u1, nu1);
267      SetCoeff(U, d, u1);
268      SetCoeff(U, 1, nu1);
269
270      GetCoeff(v0, V, 0);
271      sub(v0, v0, nv0);
272      SetCoeff(V, d-1, v0);
273      SetCoeff(V, 0, nv0);
274
275      break;
276     
277
278   case 2:
279      GetCoeff(u0, U, 0);
280      sub(u0, u0, nu0);
281      SetCoeff(U, d, u0);
282      SetCoeff(U, 0, nu0);
283
284      break;
285
286   }
287}
288
289
290void mul(zz_pXMatrix& A, zz_pXMatrix& B, zz_pXMatrix& C)
291// A = B*C, B and C are destroyed
292{
293   long db = deg(B(1,1));
294   long dc = deg(C(1,1));
295   long da = db + dc;
296
297   long k = NextPowerOfTwo(da+1);
298
299   fftRep B00, B01, B10, B11, C0, C1, T1, T2;
300   
301   TofftRep(B00, B(0,0), k); B(0,0).kill();
302   TofftRep(B01, B(0,1), k); B(0,1).kill();
303   TofftRep(B10, B(1,0), k); B(1,0).kill();
304   TofftRep(B11, B(1,1), k); B(1,1).kill();
305
306   TofftRep(C0, C(0,0), k);  C(0,0).kill();
307   TofftRep(C1, C(1,0), k);  C(1,0).kill();
308
309   mul(T1, B00, C0);
310   mul(T2, B01, C1);
311   add(T1, T1, T2);
312   FromfftRep(A(0,0), T1, 0, da);
313
314   mul(T1, B10, C0);
315   mul(T2, B11, C1);
316   add(T1, T1, T2);
317   FromfftRep(A(1,0), T1, 0, da);
318
319   TofftRep(C0, C(0,1), k);  C(0,1).kill();
320   TofftRep(C1, C(1,1), k);  C(1,1).kill();
321
322   mul(T1, B00, C0);
323   mul(T2, B01, C1);
324   add(T1, T1, T2);
325   FromfftRep(A(0,1), T1, 0, da);
326
327   mul(T1, B10, C0);
328   mul(T2, B11, C1);
329   add(T1, T1, T2);
330   FromfftRep(A(1,1), T1, 0, da);
331}
332
333void IterHalfGCD(zz_pXMatrix& M_out, zz_pX& U, zz_pX& V, long d_red)
334{
335   M_out(0,0).SetMaxLength(d_red);
336   M_out(0,1).SetMaxLength(d_red);
337   M_out(1,0).SetMaxLength(d_red);
338   M_out(1,1).SetMaxLength(d_red);
339
340   set(M_out(0,0));   clear(M_out(0,1));
341   clear(M_out(1,0)); set(M_out(1,1));
342
343   long goal = deg(U) - d_red;
344
345   if (deg(V) <= goal)
346      return;
347
348   zz_pX Q, t(INIT_SIZE, d_red);
349
350   while (deg(V) > goal) {
351      PlainDivRem(Q, U, U, V);
352      swap(U, V);
353
354      mul(t, Q, M_out(1,0));
355      sub(t, M_out(0,0), t);
356      M_out(0,0) = M_out(1,0);
357      M_out(1,0) = t;
358
359      mul(t, Q, M_out(1,1));
360      sub(t, M_out(0,1), t);
361      M_out(0,1) = M_out(1,1);
362      M_out(1,1) = t;
363   }
364}
365   
366
367
368void HalfGCD(zz_pXMatrix& M_out, const zz_pX& U, const zz_pX& V, long d_red)
369{
370   if (IsZero(V) || deg(V) <= deg(U) - d_red) {
371      set(M_out(0,0));   clear(M_out(0,1));
372      clear(M_out(1,0)); set(M_out(1,1));
373 
374      return;
375   }
376
377
378   long n = deg(U) - 2*d_red + 2;
379   if (n < 0) n = 0;
380
381   zz_pX U1, V1;
382
383   RightShift(U1, U, n);
384   RightShift(V1, V, n);
385
386   if (d_red <= NTL_zz_pX_HalfGCD_CROSSOVER) {
387      IterHalfGCD(M_out, U1, V1, d_red);
388      return;
389   }
390
391   long d1 = (d_red + 1)/2;
392   if (d1 < 1) d1 = 1;
393   if (d1 >= d_red) d1 = d_red - 1;
394
395   zz_pXMatrix M1;
396
397   HalfGCD(M1, U1, V1, d1);
398   mul(U1, V1, M1);
399
400   long d2 = deg(V1) - deg(U) + n + d_red;
401
402   if (IsZero(V1) || d2 <= 0) {
403      M_out = M1;
404      return;
405   }
406
407
408   zz_pX Q;
409   zz_pXMatrix M2;
410
411   DivRem(Q, U1, U1, V1);
412   swap(U1, V1);
413
414   HalfGCD(M2, U1, V1, d2);
415
416   zz_pX t(INIT_SIZE, deg(M1(1,1))+deg(Q)+1);
417
418   mul(t, Q, M1(1,0));
419   sub(t, M1(0,0), t);
420   swap(M1(0,0), M1(1,0));
421   swap(M1(1,0), t);
422
423   t.kill();
424
425   t.SetMaxLength(deg(M1(1,1))+deg(Q)+1);
426
427   mul(t, Q, M1(1,1));
428   sub(t, M1(0,1), t);
429   swap(M1(0,1), M1(1,1));
430   swap(M1(1,1), t);
431
432   t.kill();
433
434   mul(M_out, M2, M1); 
435}
436
437
438
439
440void XHalfGCD(zz_pXMatrix& M_out, zz_pX& U, zz_pX& V, long d_red)
441{
442   if (IsZero(V) || deg(V) <= deg(U) - d_red) {
443      set(M_out(0,0));   clear(M_out(0,1));
444      clear(M_out(1,0)); set(M_out(1,1));
445 
446      return;
447   }
448
449   long du = deg(U);
450
451   if (d_red <= NTL_zz_pX_HalfGCD_CROSSOVER) {
452      IterHalfGCD(M_out, U, V, d_red);
453      return;
454   }
455
456   long d1 = (d_red + 1)/2;
457   if (d1 < 1) d1 = 1;
458   if (d1 >= d_red) d1 = d_red - 1;
459
460   zz_pXMatrix M1;
461
462   HalfGCD(M1, U, V, d1);
463   mul(U, V, M1);
464
465   long d2 = deg(V) - du + d_red;
466
467   if (IsZero(V) || d2 <= 0) {
468      M_out = M1;
469      return;
470   }
471
472
473   zz_pX Q;
474   zz_pXMatrix M2;
475
476   DivRem(Q, U, U, V);
477   swap(U, V);
478
479   XHalfGCD(M2, U, V, d2);
480
481   zz_pX t(INIT_SIZE, deg(M1(1,1))+deg(Q)+1);
482
483   mul(t, Q, M1(1,0));
484   sub(t, M1(0,0), t);
485   swap(M1(0,0), M1(1,0));
486   swap(M1(1,0), t);
487
488   t.kill();
489
490   t.SetMaxLength(deg(M1(1,1))+deg(Q)+1);
491
492   mul(t, Q, M1(1,1));
493   sub(t, M1(0,1), t);
494   swap(M1(0,1), M1(1,1));
495   swap(M1(1,1), t);
496
497   t.kill();
498
499   mul(M_out, M2, M1); 
500}
501
502void HalfGCD(zz_pX& U, zz_pX& V)
503{
504   long d_red = (deg(U)+1)/2;
505
506   if (IsZero(V) || deg(V) <= deg(U) - d_red) {
507      return;
508   }
509
510   long du = deg(U);
511
512
513   long d1 = (d_red + 1)/2;
514   if (d1 < 1) d1 = 1;
515   if (d1 >= d_red) d1 = d_red - 1;
516
517   zz_pXMatrix M1;
518
519   HalfGCD(M1, U, V, d1);
520   mul(U, V, M1);
521
522   long d2 = deg(V) - du + d_red;
523
524   if (IsZero(V) || d2 <= 0) {
525      return;
526   }
527
528   M1(0,0).kill();
529   M1(0,1).kill();
530   M1(1,0).kill();
531   M1(1,1).kill();
532
533
534   zz_pX Q;
535
536   DivRem(Q, U, U, V);
537   swap(U, V);
538
539   HalfGCD(M1, U, V, d2);
540
541   mul(U, V, M1); 
542}
543
544
545void GCD(zz_pX& d, const zz_pX& u, const zz_pX& v)
546{
547   zz_pX u1, v1;
548
549   u1 = u;
550   v1 = v;
551
552   if (deg(u1) == deg(v1)) {
553      if (IsZero(u1)) {
554         clear(d);
555         return;
556      }
557
558      rem(v1, v1, u1);
559   }
560   else if (deg(u1) < deg(v1)) {
561      swap(u1, v1);
562   }
563
564   // deg(u1) > deg(v1)
565
566   while (deg(u1) > NTL_zz_pX_GCD_CROSSOVER && !IsZero(v1)) {
567      HalfGCD(u1, v1);
568
569      if (!IsZero(v1)) {
570         rem(u1, u1, v1);
571         swap(u1, v1);
572      }
573   }
574
575   PlainGCD(d, u1, v1);
576}
577
578
579
580void XGCD(zz_pX& d, zz_pX& s, zz_pX& t, const zz_pX& a, const zz_pX& b)
581{
582   zz_p w;
583
584   if (IsZero(a) && IsZero(b)) {
585      clear(d);
586      set(s);
587      clear(t);
588      return;
589   }
590
591   zz_pX U, V, Q;
592
593   U = a;
594   V = b;
595
596   long flag = 0;
597
598   if (deg(U) == deg(V)) {
599      DivRem(Q, U, U, V);
600      swap(U, V);
601      flag = 1;
602   }
603   else if (deg(U) < deg(V)) {
604      swap(U, V);
605      flag = 2;
606   }
607
608   zz_pXMatrix M;
609
610   XHalfGCD(M, U, V, deg(U)+1);
611
612   d = U;
613
614   if (flag == 0) {
615      s = M(0,0); 
616      t = M(0,1);
617   }
618   else if (flag == 1) {
619      s = M(0,1);
620      mul(t, Q, M(0,1));
621      sub(t, M(0,0), t);
622   }
623   else {  /* flag == 2 */
624      s = M(0,1);
625      t = M(0,0);
626   }
627
628   // normalize
629
630   inv(w, LeadCoeff(d));
631   mul(d, d, w);
632   mul(s, s, w);
633   mul(t, t, w);
634}
635
636     
637
638
639
640
641
642void IterBuild(zz_p* a, long n)
643{
644   long i, k;
645   zz_p b, t;
646
647   if (n <= 0) return;
648
649   negate(a[0], a[0]);
650
651   for (k = 1; k <= n-1; k++) {
652      negate(b, a[k]);
653      add(a[k], b, a[k-1]);
654      for (i = k-1; i >= 1; i--) {
655         mul(t, a[i], b);
656         add(a[i], t, a[i-1]);
657      }
658      mul(a[0], a[0], b);
659   }
660} 
661
662void mul(zz_p* x, const zz_p* a, const zz_p* b, long n)
663{
664   zz_p t, accum;
665
666   long i, j, jmin, jmax;
667
668   long d = 2*n-1;
669
670   for (i = 0; i <= d; i++) {
671      jmin = max(0, i-(n-1));
672      jmax = min(n-1, i);
673      clear(accum);
674      for (j = jmin; j <= jmax; j++) {
675         mul(t, (a[j]), (b[i-j]));
676         add(accum, accum, t);
677      }
678      if (i >= n) {
679         add(accum, accum, (a[i-n]));
680         add(accum, accum, (b[i-n]));
681      }
682
683      x[i] = accum;
684   }
685}
686
687
688void BuildFromRoots(zz_pX& x, const vec_zz_p& a)
689{
690   long n = a.length();
691
692   if (n == 0) {
693      set(x);
694      return;
695   }
696
697   long k0 = NextPowerOfTwo(NTL_zz_pX_MUL_CROSSOVER)-1;
698   long crossover = 1L << k0;
699
700   if (n <= NTL_zz_pX_MUL_CROSSOVER) {
701      x.rep.SetMaxLength(n+1);
702      x.rep = a;
703      IterBuild(&x.rep[0], n);
704      x.rep.SetLength(n+1);
705      SetCoeff(x, n);
706      return;
707   }
708
709   long k = NextPowerOfTwo(n);
710
711   long m = 1L << k;
712   long i, j;
713   long l, width;
714
715   zz_pX b(INIT_SIZE, m+1);
716
717   b.rep = a;
718   b.rep.SetLength(m+1);
719   for (i = n; i < m; i++)
720      clear(b.rep[i]);
721
722   set(b.rep[m]);
723   
724   fftRep R1(INIT_SIZE, k), R2(INIT_SIZE, k);
725
726
727   zz_p t1, one;
728   set(one);
729
730   vec_zz_p G(INIT_SIZE, crossover), H(INIT_SIZE, crossover);
731   zz_p *g = G.elts();
732   zz_p *h = H.elts();
733   zz_p *tmp;
734   
735   for (i = 0; i < m; i+= crossover) {
736      for (j = 0; j < crossover; j++)
737         negate(g[j], b.rep[i+j]);
738
739      if (k0 > 0) {
740         for (j = 0; j < crossover; j+=2) {
741            mul(t1, g[j], g[j+1]);
742            add(g[j+1], g[j], g[j+1]);
743            g[j] = t1;
744         }
745      }
746   
747      for (l = 1; l < k0; l++) {
748         width = 1L << l;
749
750         for (j = 0; j < crossover; j += 2*width)
751            mul(&h[j], &g[j], &g[j+width], width);
752     
753         tmp = g; g = h; h = tmp;
754      }
755
756      for (j = 0; j < crossover; j++)
757         b.rep[i+j] = g[j];
758   }
759
760   for (l = k0; l < k; l++) {
761      width = 1L << l;
762      for (i = 0; i < m; i += 2*width) {
763         t1 = b.rep[i+width];
764         set(b.rep[i+width]);
765         TofftRep(R1, b, l+1, i, i+width);
766         b.rep[i+width] = t1;
767         t1 = b.rep[i+2*width];
768         set(b.rep[i+2*width]);
769         TofftRep(R2, b, l+1, i+width, i+2*width);
770         b.rep[i+2*width] = t1;
771         mul(R1, R1, R2);
772         FromfftRep(&b.rep[i], R1, 0, 2*width-1);
773         sub(b.rep[i], b.rep[i], one);
774      }
775   }
776
777   x.rep.SetLength(n+1);
778   long delta = m-n;
779   for (i = 0; i <= n; i++)
780     x.rep[i] = b.rep[i+delta];
781
782   // no need to normalize
783}
784
785
786
787void eval(zz_p& b, const zz_pX& f, zz_p a)
788// does a Horner evaluation
789{
790   zz_p acc;
791   long i;
792
793   clear(acc);
794   for (i = deg(f); i >= 0; i--) {
795      mul(acc, acc, a);
796      add(acc, acc, f.rep[i]);
797   }
798
799   b = acc;
800}
801
802
803
804void eval(vec_zz_p& b, const zz_pX& f, const vec_zz_p& a)
805// naive algorithm:  repeats Horner
806{
807   if (&b == &f.rep) {
808      vec_zz_p bb;
809      eval(bb, f, a);
810      b = bb;
811      return;
812   }
813
814   long m = a.length();
815   b.SetLength(m);
816   long i;
817   for (i = 0; i < m; i++) 
818      eval(b[i], f, a[i]);
819}
820
821
822
823
824void interpolate(zz_pX& f, const vec_zz_p& a, const vec_zz_p& b)
825{
826   long m = a.length();
827   if (b.length() != m) Error("interpolate: vector length mismatch");
828
829   if (m == 0) {
830      clear(f);
831      return;
832   }
833
834   vec_zz_p prod;
835   prod = a;
836
837   zz_p t1, t2;
838
839   long k, i;
840
841   vec_zz_p res;
842   res.SetLength(m);
843
844   for (k = 0; k < m; k++) {
845
846      const zz_p& aa = a[k];
847
848      set(t1);
849      for (i = k-1; i >= 0; i--) {
850         mul(t1, t1, aa);
851         add(t1, t1, prod[i]);
852      }
853
854      clear(t2);
855      for (i = k-1; i >= 0; i--) {
856         mul(t2, t2, aa);
857         add(t2, t2, res[i]);
858      }
859
860
861      inv(t1, t1);
862      sub(t2, b[k], t2);
863      mul(t1, t1, t2);
864
865      for (i = 0; i < k; i++) {
866         mul(t2, prod[i], t1);
867         add(res[i], res[i], t2);
868      }
869
870      res[k] = t1;
871
872      if (k < m-1) {
873         if (k == 0)
874            negate(prod[0], prod[0]);
875         else {
876            negate(t1, a[k]);
877            add(prod[k], t1, prod[k-1]);
878            for (i = k-1; i >= 1; i--) {
879               mul(t2, prod[i], t1);
880               add(prod[i], t2, prod[i-1]);
881            }
882            mul(prod[0], prod[0], t1);
883         }
884      }
885   }
886
887   while (m > 0 && IsZero(res[m-1])) m--;
888   res.SetLength(m);
889   f.rep = res;
890}
891
892NTL_vector_impl(zz_pX,vec_zz_pX)
893
894NTL_eq_vector_impl(zz_pX,vec_zz_pX)
895
896
897
898   
899void InnerProduct(zz_pX& x, const vec_zz_p& v, long low, long high, 
900                   const vec_zz_pX& H, long n, vec_zz_p& t)
901{
902   zz_p s;
903   long i, j;
904
905   zz_p *tp = t.elts();
906
907   for (j = 0; j < n; j++)
908      clear(tp[j]);
909
910
911   long p = zz_p::modulus();
912   double pinv = zz_p::ModulusInverse();
913
914   high = min(high, v.length()-1);
915   for (i = low; i <= high; i++) {
916      const vec_zz_p& h = H[i-low].rep;
917      long m = h.length();
918      zz_p w = (v[i]);
919
920      long W = rep(w);
921      mulmod_precon_t Wpinv = PrepMulModPrecon(W, p, pinv); // ((double) W)*pinv;
922      const zz_p *hp = h.elts();
923
924      for (j = 0; j < m; j++) {
925         long S = MulModPrecon(rep(hp[j]), W, p, Wpinv);
926         S = AddMod(S, rep(tp[j]), p);
927         tp[j].LoopHole() = S;
928      }
929   }
930
931   x.rep = t;
932   x.normalize();
933}
934
935
936void CompMod(zz_pX& x, const zz_pX& g, const zz_pXArgument& A, 
937             const zz_pXModulus& F)
938{
939   if (deg(g) <= 0) {
940      x = g;
941      return;
942   }
943
944
945   zz_pX s, t;
946   vec_zz_p scratch(INIT_SIZE, F.n);
947
948   long m = A.H.length() - 1;
949   long l = ((g.rep.length()+m-1)/m) - 1;
950
951   zz_pXMultiplier M;
952   build(M, A.H[m], F);
953
954   InnerProduct(t, g.rep, l*m, l*m + m - 1, A.H, F.n, scratch);
955   for (long i = l-1; i >= 0; i--) {
956      InnerProduct(s, g.rep, i*m, i*m + m - 1, A.H, F.n, scratch);
957      MulMod(t, t, M, F);
958      add(t, t, s);
959   }
960
961   x = t;
962}
963
964
965void build(zz_pXArgument& A, const zz_pX& h, const zz_pXModulus& F, long m)
966{
967   if (m <= 0 || deg(h) >= F.n) Error("build: bad args");
968
969   if (m > F.n) m = F.n;
970
971   long i;
972
973   if (zz_pXArgBound > 0) {
974      double sz = 1;
975      sz = sz*F.n;
976      sz = sz+6;
977      sz = sz*(sizeof (long));
978      sz = sz/1024;
979      m = min(m, long(zz_pXArgBound/sz));
980      m = max(m, 1);
981   }
982
983   zz_pXMultiplier M;
984
985   build(M, h, F);
986
987   A.H.SetLength(m+1);
988
989   set(A.H[0]);
990   A.H[1] = h;
991   for (i = 2; i <= m; i++) 
992      MulMod(A.H[i], A.H[i-1], M, F);
993}
994
995
996
997
998long zz_pXArgBound = 0;
999
1000
1001void CompMod(zz_pX& x, const zz_pX& g, const zz_pX& h, const zz_pXModulus& F)
1002   // x = g(h) mod f
1003{
1004   long m = SqrRoot(g.rep.length());
1005
1006   if (m == 0) {
1007      clear(x);
1008      return;
1009   }
1010
1011   zz_pXArgument A;
1012
1013   build(A, h, F, m);
1014
1015   CompMod(x, g, A, F);
1016}
1017
1018
1019
1020
1021void Comp2Mod(zz_pX& x1, zz_pX& x2, const zz_pX& g1, const zz_pX& g2,
1022              const zz_pX& h, const zz_pXModulus& F)
1023
1024{
1025   long m = SqrRoot(g1.rep.length() + g2.rep.length());
1026
1027   if (m == 0) {
1028      clear(x1);
1029      clear(x2);
1030      return;
1031   }
1032
1033   zz_pXArgument A;
1034
1035   build(A, h, F, m);
1036
1037   zz_pX xx1, xx2;
1038
1039   CompMod(xx1, g1, A, F);
1040   CompMod(xx2, g2, A, F);
1041
1042   x1 = xx1;
1043   x2 = xx2;
1044}
1045
1046void Comp3Mod(zz_pX& x1, zz_pX& x2, zz_pX& x3, 
1047              const zz_pX& g1, const zz_pX& g2, const zz_pX& g3,
1048              const zz_pX& h, const zz_pXModulus& F)
1049
1050{
1051   long m = SqrRoot(g1.rep.length() + g2.rep.length() + g3.rep.length());
1052
1053   if (m == 0) {
1054      clear(x1);
1055      clear(x2);
1056      clear(x3);
1057      return;
1058   }
1059
1060   zz_pXArgument A;
1061
1062   build(A, h, F, m);
1063
1064   zz_pX xx1, xx2, xx3;
1065
1066   CompMod(xx1, g1, A, F);
1067   CompMod(xx2, g2, A, F);
1068   CompMod(xx3, g3, A, F);
1069
1070   x1 = xx1;
1071   x2 = xx2;
1072   x3 = xx3;
1073}
1074
1075static void StripZeroes(vec_zz_p& x)
1076{
1077   long n = x.length();
1078   while (n > 0 && IsZero(x[n-1]))
1079      n--;
1080   x.SetLength(n);
1081}
1082
1083
1084void PlainUpdateMap(vec_zz_p& xx, const vec_zz_p& a, 
1085                    const zz_pX& b, const zz_pX& f)
1086{
1087   long n = deg(f);
1088   long i, m;
1089
1090   if (IsZero(b)) {
1091      xx.SetLength(0);
1092      return;
1093   }
1094
1095   m = n-1 - deg(b);
1096
1097   vec_zz_p x(INIT_SIZE, n);
1098
1099   for (i = 0; i <= m; i++)
1100      InnerProduct(x[i], a, b.rep, i);
1101
1102   if (deg(b) != 0) {
1103      zz_pX c(INIT_SIZE, n);
1104      LeftShift(c, b, m);
1105
1106      for (i = m+1; i < n; i++) {
1107         MulByXMod(c, c, f);
1108         InnerProduct(x[i], a, c.rep);
1109      }
1110   }
1111
1112   xx = x;
1113}
1114   
1115
1116
1117
1118void UpdateMap(vec_zz_p& x, const vec_zz_p& aa, 
1119               const zz_pXMultiplier& B, const zz_pXModulus& F)
1120{
1121   long n = F.n;
1122
1123   vec_zz_p a;
1124   a = aa;
1125   StripZeroes(a);
1126
1127   if (a.length() > n) Error("UpdateMap: bad args");
1128   long i;
1129
1130   if (!B.UseFFT) {
1131      PlainUpdateMap(x, a, B.b, F.f);
1132      StripZeroes(x);
1133      return;
1134   }
1135
1136   fftRep R1(INIT_SIZE, F.k), R2(INIT_SIZE, F.l);
1137   vec_zz_p V1(INIT_SIZE, n);
1138
1139
1140   RevTofftRep(R1, a, F.k, 0, a.length()-1, 0);
1141   mul(R2, R1, F.FRep);
1142   RevFromfftRep(V1, R2, 0, n-2);
1143   for (i = 0; i <= n-2; i++)  negate(V1[i], V1[i]);
1144   RevTofftRep(R2, V1, F.l, 0, n-2, n-1);
1145   mul(R2, R2, B.B1);
1146   mul(R1, R1, B.B2);
1147
1148   AddExpand(R2, R1);
1149   RevFromfftRep(x, R2, 0, n-1);
1150   StripZeroes(x);
1151}
1152
1153   
1154
1155void ProjectPowers(vec_zz_p& x, const vec_zz_p& a, long k,
1156                   const zz_pXArgument& H, const zz_pXModulus& F)
1157
1158{
1159   long n = F.n;
1160
1161   if (a.length() > n || k < 0 || NTL_OVERFLOW(k, 1, 0))
1162      Error("ProjectPowers: bad args");
1163
1164   long m = H.H.length()-1;
1165   long l = (k+m-1)/m - 1;
1166
1167   zz_pXMultiplier M;
1168   build(M, H.H[m], F);
1169
1170   vec_zz_p s(INIT_SIZE, n);
1171   s = a;
1172   StripZeroes(s);
1173
1174   x.SetLength(k);
1175
1176   for (long i = 0; i <= l; i++) {
1177      long m1 = min(m, k-i*m);
1178      zz_p* w = &x[i*m];
1179      for (long j = 0; j < m1; j++)
1180         InnerProduct(w[j], H.H[j].rep, s);
1181      if (i < l)
1182         UpdateMap(s, s, M, F);
1183   }
1184}
1185
1186
1187
1188void ProjectPowers(vec_zz_p& x, const vec_zz_p& a, long k,
1189                   const zz_pX& h, const zz_pXModulus& F)
1190
1191{
1192   if (a.length() > F.n || k < 0) Error("ProjectPowers: bad args");
1193
1194   if (k == 0) {
1195      x.SetLength(0);
1196      return;
1197   }
1198
1199   long m = SqrRoot(k);
1200
1201   zz_pXArgument H;
1202
1203   build(H, h, F, m);
1204   ProjectPowers(x, a, k, H, F);
1205}
1206
1207
1208void BerlekampMassey(zz_pX& h, const vec_zz_p& a, long m)
1209{
1210   zz_pX Lambda, Sigma, Temp;
1211   long L;
1212   zz_p Delta, Delta1, t1;
1213   long shamt;
1214
1215   // cerr << "*** " << m << "\n";
1216
1217   Lambda.SetMaxLength(m+1);
1218   Sigma.SetMaxLength(m+1);
1219   Temp.SetMaxLength(m+1);
1220
1221   L = 0;
1222   set(Lambda);
1223   clear(Sigma);
1224   set(Delta);
1225   shamt = 0;
1226
1227   long i, r, dl;
1228
1229   for (r = 1; r <= 2*m; r++) {
1230      // cerr << r << "--";
1231      clear(Delta1);
1232      dl = deg(Lambda);
1233      for (i = 0; i <= dl; i++) {
1234         mul(t1, Lambda.rep[i], a[r-i-1]);
1235         add(Delta1, Delta1, t1);
1236      }
1237
1238      if (IsZero(Delta1)) {
1239         shamt++;
1240         // cerr << "case 1: " << deg(Lambda) << " " << deg(Sigma) << " " << shamt << "\n";
1241      }
1242      else if (2*L < r) {
1243         div(t1, Delta1, Delta);
1244         mul(Temp, Sigma, t1);
1245         Sigma = Lambda;
1246         ShiftSub(Lambda, Temp, shamt+1);
1247         shamt = 0;
1248         L = r-L;
1249         Delta = Delta1;
1250         // cerr << "case 2: " << deg(Lambda) << " " << deg(Sigma) << " " << shamt << "\n";
1251      }
1252      else {
1253         shamt++;
1254         div(t1, Delta1, Delta);
1255         mul(Temp, Sigma, t1);
1256         ShiftSub(Lambda, Temp, shamt);
1257         // cerr << "case 3: " << deg(Lambda) << " " << deg(Sigma) << " " << shamt << "\n";
1258      }
1259   }
1260
1261   // cerr << "finished: " << L << " " << deg(Lambda) << "\n";
1262
1263   dl = deg(Lambda);
1264   h.rep.SetLength(L + 1);
1265
1266   for (i = 0; i < L - dl; i++)
1267      clear(h.rep[i]);
1268
1269   for (i = L - dl; i <= L; i++)
1270      h.rep[i] = Lambda.rep[L - i];
1271}
1272
1273
1274void GCDMinPolySeq(zz_pX& h, const vec_zz_p& x, long m)
1275{
1276   long i;
1277   zz_pX a, b;
1278   zz_pXMatrix M;
1279   zz_p t;
1280
1281   a.rep.SetLength(2*m);
1282   for (i = 0; i < 2*m; i++) a.rep[i] = x[2*m-1-i];
1283   a.normalize();
1284
1285   SetCoeff(b, 2*m);
1286
1287   HalfGCD(M, b, a, m+1);
1288
1289   /* make monic */
1290
1291   inv(t, LeadCoeff(M(1,1)));
1292   mul(h, M(1,1), t);
1293}
1294
1295
1296void MinPolySeq(zz_pX& h, const vec_zz_p& a, long m)
1297{
1298   if (m < 0 || NTL_OVERFLOW(m, 1, 0)) Error("MinPoly: bad args");
1299   if (a.length() < 2*m) Error("MinPoly: sequence too short");
1300
1301   if (m > NTL_zz_pX_BERMASS_CROSSOVER)
1302      GCDMinPolySeq(h, a, m);
1303   else
1304      BerlekampMassey(h, a, m);
1305}
1306
1307
1308void DoMinPolyMod(zz_pX& h, const zz_pX& g, const zz_pXModulus& F, long m,
1309               const vec_zz_p& R) 
1310{
1311   vec_zz_p x;
1312
1313   ProjectPowers(x, R, 2*m, g, F);
1314   MinPolySeq(h, x, m);
1315}
1316
1317
1318void ProbMinPolyMod(zz_pX& h, const zz_pX& g, const zz_pXModulus& F, long m)
1319{
1320   long n = F.n;
1321   if (m < 1 || m > n) Error("ProbMinPoly: bad args");
1322
1323   long i;
1324   vec_zz_p R(INIT_SIZE, n);
1325
1326   for (i = 0; i < n; i++) random(R[i]);
1327   DoMinPolyMod(h, g, F, m, R);
1328}
1329
1330void MinPolyMod(zz_pX& hh, const zz_pX& g, const zz_pXModulus& F, long m)
1331{
1332   zz_pX h, h1;
1333   long n = F.n;
1334   if (m < 1 || m > n) Error("MinPoly: bad args");
1335
1336   /* probabilistically compute min-poly */
1337
1338   ProbMinPolyMod(h, g, F, m);
1339   if (deg(h) == m) { hh = h; return; }
1340   CompMod(h1, h, g, F);
1341   if (IsZero(h1)) { hh = h; return; }
1342
1343   /* not completely successful...must iterate */
1344
1345   long i;
1346
1347   zz_pX h2, h3;
1348   zz_pXMultiplier H1;
1349   vec_zz_p R(INIT_SIZE, n);
1350
1351   for (;;) {
1352      R.SetLength(n);
1353      for (i = 0; i < n; i++) random(R[i]);
1354      build(H1, h1, F);
1355      UpdateMap(R, R, H1, F);
1356      DoMinPolyMod(h2, g, F, m-deg(h), R);
1357
1358      mul(h, h, h2);
1359      if (deg(h) == m) { hh = h; return; }
1360      CompMod(h3, h2, g, F);
1361      MulMod(h1, h3, H1, F);
1362      if (IsZero(h1)) { hh = h; return; }
1363   }
1364}
1365
1366void IrredPolyMod(zz_pX& h, const zz_pX& g, const zz_pXModulus& F, long m)
1367{
1368   vec_zz_p R(INIT_SIZE, 1);
1369   if (m < 1 || m > F.n) Error("IrredPoly: bad args");
1370
1371   set(R[0]);
1372   DoMinPolyMod(h, g, F, m, R);
1373}
1374
1375
1376
1377void diff(zz_pX& x, const zz_pX& a)
1378{
1379   long n = deg(a);
1380   long i;
1381
1382   if (n <= 0) {
1383      clear(x);
1384      return;
1385   }
1386
1387   if (&x != &a)
1388      x.rep.SetLength(n);
1389
1390   for (i = 0; i <= n-1; i++) {
1391      mul(x.rep[i], a.rep[i+1], i+1);
1392   }
1393
1394   if (&x == &a)
1395      x.rep.SetLength(n);
1396
1397   x.normalize();
1398}
1399
1400void MakeMonic(zz_pX& x)
1401{
1402   if (IsZero(x))
1403      return;
1404
1405   if (IsOne(LeadCoeff(x)))
1406      return;
1407
1408   zz_p t;
1409
1410   inv(t, LeadCoeff(x));
1411   mul(x, x, t);
1412}
1413
1414
1415
1416
1417     
1418void PlainMulTrunc(zz_pX& x, const zz_pX& a, const zz_pX& b, long n)
1419{
1420   zz_pX y;
1421   mul(y, a, b);
1422   trunc(x, y, n);
1423}
1424
1425
1426void FFTMulTrunc(zz_pX& x, const zz_pX& a, const zz_pX& b, long n)
1427{
1428   if (IsZero(a) || IsZero(b)) {
1429      clear(x);
1430      return;
1431   }
1432
1433   long d = deg(a) + deg(b);
1434   if (n > d + 1)
1435      n = d + 1;
1436
1437   long k = NextPowerOfTwo(d + 1);
1438   fftRep R1(INIT_SIZE, k), R2(INIT_SIZE, k);
1439
1440   TofftRep(R1, a, k);
1441   TofftRep(R2, b, k);
1442   mul(R1, R1, R2);
1443   FromfftRep(x, R1, 0, n-1);
1444}
1445
1446void MulTrunc(zz_pX& x, const zz_pX& a, const zz_pX& b, long n)
1447{
1448   if (n < 0) Error("MulTrunc: bad args");
1449
1450   if (deg(a) <= NTL_zz_pX_MUL_CROSSOVER || deg(b) <= NTL_zz_pX_MUL_CROSSOVER)
1451      PlainMulTrunc(x, a, b, n);
1452   else
1453      FFTMulTrunc(x, a, b, n);
1454}
1455
1456void PlainSqrTrunc(zz_pX& x, const zz_pX& a, long n)
1457{
1458   zz_pX y;
1459   sqr(y, a);
1460   trunc(x, y, n);
1461}
1462
1463
1464void FFTSqrTrunc(zz_pX& x, const zz_pX& a, long n)
1465{
1466   if (IsZero(a)) {
1467      clear(x);
1468      return;
1469   }
1470
1471   long d = 2*deg(a);
1472   if (n > d + 1)
1473      n = d + 1;
1474
1475   long k = NextPowerOfTwo(d + 1);
1476   fftRep R1(INIT_SIZE, k);
1477
1478   TofftRep(R1, a, k);
1479   mul(R1, R1, R1);
1480   FromfftRep(x, R1, 0, n-1);
1481}
1482
1483void SqrTrunc(zz_pX& x, const zz_pX& a, long n)
1484{
1485   if (n < 0) Error("SqrTrunc: bad args");
1486
1487   if (deg(a) <= NTL_zz_pX_MUL_CROSSOVER)
1488      PlainSqrTrunc(x, a, n);
1489   else
1490      FFTSqrTrunc(x, a, n);
1491}
1492
1493
1494
1495void FastTraceVec(vec_zz_p& S, const zz_pX& f)
1496{
1497   long n = deg(f);
1498
1499   if (n <= 0) 
1500      Error("FastTraceVec: bad args");
1501
1502   if (n == 0) {
1503      S.SetLength(0);
1504      return;
1505   }
1506
1507   if (n == 1) {
1508      S.SetLength(1);
1509      set(S[0]);
1510      return;
1511   }
1512   
1513   long i;
1514   zz_pX f1;
1515
1516   f1.rep.SetLength(n-1);
1517   for (i = 0; i <= n-2; i++)
1518      f1.rep[i] = f.rep[n-i];
1519   f1.normalize();
1520
1521   zz_pX f2;
1522   f2.rep.SetLength(n-1);
1523   for (i = 0; i <= n-2; i++)
1524      mul(f2.rep[i], f.rep[n-1-i], i+1);
1525   f2.normalize();
1526
1527   zz_pX f3;
1528   InvTrunc(f3, f1, n-1);
1529   MulTrunc(f3, f3, f2, n-1);
1530
1531   S.SetLength(n);
1532
1533   S[0] = n;
1534   for (i = 1; i < n; i++)
1535      negate(S[i], coeff(f3, i-1));
1536}
1537
1538
1539void PlainTraceVec(vec_zz_p& S, const zz_pX& ff)
1540{
1541   if (deg(ff) <= 0)
1542      Error("TraceVec: bad args");
1543
1544   zz_pX f;
1545   f = ff;
1546
1547   MakeMonic(f);
1548
1549   long n = deg(f);
1550
1551   S.SetLength(n);
1552
1553   if (n == 0)
1554      return;
1555
1556   long k, i;
1557   zz_p acc, t;
1558
1559   const zz_p *fp = f.rep.elts();;
1560   zz_p *sp = S.elts();
1561
1562   sp[0] = n;
1563
1564   for (k = 1; k < n; k++) {
1565      mul(acc, fp[n-k], k);
1566
1567      for (i = 1; i < k; i++) {
1568         mul(t, fp[n-i], rep(sp[k-i]));
1569         add(acc, acc, t);
1570      }
1571
1572      negate(sp[k], acc);
1573   }
1574}
1575
1576void TraceVec(vec_zz_p& S, const zz_pX& f)
1577{
1578   if (deg(f) <= NTL_zz_pX_TRACE_CROSSOVER)
1579      PlainTraceVec(S, f);
1580   else
1581      FastTraceVec(S, f);
1582}
1583
1584void ComputeTraceVec(const zz_pXModulus& F)
1585{
1586   vec_zz_p& S = *((vec_zz_p *) &F.tracevec);
1587
1588   if (S.length() > 0)
1589      return;
1590
1591   if (!F.UseFFT) {
1592      PlainTraceVec(S, F.f);
1593      return;
1594   }
1595
1596   long i;
1597   long n = F.n;
1598
1599   fftRep R;
1600   zz_pX P, g;
1601
1602   g.rep.SetLength(n-1);
1603   for (i = 1; i < n; i++)
1604      mul(g.rep[n-i-1], F.f.rep[n-i], i); 
1605   g.normalize();
1606
1607   TofftRep(R, g, F.l);
1608   mul(R, R, F.HRep);
1609   FromfftRep(P, R, n-2, 2*n-4);
1610
1611   S.SetLength(n);
1612
1613   S[0] = n;
1614   for (i = 1; i < n; i++)
1615      negate(S[i], coeff(P, n-1-i));
1616}
1617
1618void TraceMod(zz_p& x, const zz_pX& a, const zz_pXModulus& F)
1619{
1620   long n = F.n;
1621
1622   if (deg(a) >= n)
1623      Error("trace: bad args");
1624
1625   if (F.tracevec.length() == 0) 
1626      TraceVec(F);
1627
1628   InnerProduct(x, a.rep, F.tracevec);
1629}
1630
1631
1632void TraceMod(zz_p& x, const zz_pX& a, const zz_pX& f)
1633{
1634   if (deg(a) >= deg(f) || deg(f) <= 0)
1635      Error("trace: bad args");
1636
1637   project(x, TraceVec(f), a);
1638}
1639
1640
1641void PlainResultant(zz_p& rres, const zz_pX& a, const zz_pX& b)
1642{
1643   zz_p res;
1644 
1645   if (IsZero(a) || IsZero(b))
1646      clear(res);
1647   else if (deg(a) == 0 && deg(b) == 0) 
1648      set(res);
1649   else {
1650      long d0, d1, d2;
1651      zz_p lc;
1652      set(res);
1653
1654      long n = max(deg(a),deg(b)) + 1;
1655      zz_pX u(INIT_SIZE, n), v(INIT_SIZE, n);
1656
1657      u = a;
1658      v = b;
1659
1660      for (;;) {
1661         d0 = deg(u);
1662         d1 = deg(v);
1663         lc = LeadCoeff(v);
1664
1665         PlainRem(u, u, v);
1666         swap(u, v);
1667
1668         d2 = deg(v);
1669         if (d2 >= 0) {
1670            power(lc, lc, d0-d2);
1671            mul(res, res, lc);
1672            if (d0 & d1 & 1) negate(res, res);
1673         }
1674         else {
1675            if (d1 == 0) {
1676               power(lc, lc, d0);
1677               mul(res, res, lc);
1678            }
1679            else
1680               clear(res);
1681       
1682            break;
1683         }
1684      }
1685
1686      rres = res;
1687   }
1688}
1689
1690
1691void ResIterHalfGCD(zz_pXMatrix& M_out, zz_pX& U, zz_pX& V, long d_red,
1692                    vec_zz_p& cvec, vec_long& dvec)
1693{
1694   M_out(0,0).SetMaxLength(d_red);
1695   M_out(0,1).SetMaxLength(d_red);
1696   M_out(1,0).SetMaxLength(d_red);
1697   M_out(1,1).SetMaxLength(d_red);
1698
1699   set(M_out(0,0));   clear(M_out(0,1));
1700   clear(M_out(1,0)); set(M_out(1,1));
1701
1702   long goal = deg(U) - d_red;
1703
1704   if (deg(V) <= goal)
1705      return;
1706
1707   zz_pX Q, t(INIT_SIZE, d_red);
1708
1709
1710   while (deg(V) > goal) {
1711      append(cvec, LeadCoeff(V));
1712      append(dvec, dvec[dvec.length()-1]-deg(U)+deg(V));
1713      PlainDivRem(Q, U, U, V);
1714      swap(U, V);
1715
1716      mul(t, Q, M_out(1,0));
1717      sub(t, M_out(0,0), t);
1718      M_out(0,0) = M_out(1,0);
1719      M_out(1,0) = t;
1720
1721      mul(t, Q, M_out(1,1));
1722      sub(t, M_out(0,1), t);
1723      M_out(0,1) = M_out(1,1);
1724      M_out(1,1) = t;
1725   }
1726}
1727   
1728
1729
1730void ResHalfGCD(zz_pXMatrix& M_out, const zz_pX& U, const zz_pX& V, long d_red,
1731                vec_zz_p& cvec, vec_long& dvec)
1732{
1733   if (IsZero(V) || deg(V) <= deg(U) - d_red) {
1734      set(M_out(0,0));   clear(M_out(0,1));
1735      clear(M_out(1,0)); set(M_out(1,1));
1736 
1737      return;
1738   }
1739
1740
1741   long n = deg(U) - 2*d_red + 2;
1742   if (n < 0) n = 0;
1743
1744   zz_pX U1, V1;
1745
1746   RightShift(U1, U, n);
1747   RightShift(V1, V, n);
1748
1749   if (d_red <= NTL_zz_pX_HalfGCD_CROSSOVER) { 
1750      ResIterHalfGCD(M_out, U1, V1, d_red, cvec, dvec);
1751      return;
1752   }
1753
1754   long d1 = (d_red + 1)/2;
1755   if (d1 < 1) d1 = 1;
1756   if (d1 >= d_red) d1 = d_red - 1;
1757
1758   zz_pXMatrix M1;
1759
1760   ResHalfGCD(M1, U1, V1, d1, cvec, dvec);
1761   mul(U1, V1, M1);
1762
1763   long d2 = deg(V1) - deg(U) + n + d_red;
1764
1765   if (IsZero(V1) || d2 <= 0) {
1766      M_out = M1;
1767      return;
1768   }
1769
1770
1771   zz_pX Q;
1772   zz_pXMatrix M2;
1773
1774   append(cvec, LeadCoeff(V1));
1775   append(dvec, dvec[dvec.length()-1]-deg(U1)+deg(V1));
1776   DivRem(Q, U1, U1, V1);
1777   swap(U1, V1);
1778
1779   ResHalfGCD(M2, U1, V1, d2, cvec, dvec);
1780
1781   zz_pX t(INIT_SIZE, deg(M1(1,1))+deg(Q)+1);
1782
1783   mul(t, Q, M1(1,0));
1784   sub(t, M1(0,0), t);
1785   swap(M1(0,0), M1(1,0));
1786   swap(M1(1,0), t);
1787
1788   t.kill();
1789
1790   t.SetMaxLength(deg(M1(1,1))+deg(Q)+1);
1791
1792   mul(t, Q, M1(1,1));
1793   sub(t, M1(0,1), t);
1794   swap(M1(0,1), M1(1,1));
1795   swap(M1(1,1), t);
1796
1797   t.kill();
1798
1799   mul(M_out, M2, M1); 
1800}
1801
1802void ResHalfGCD(zz_pX& U, zz_pX& V, vec_zz_p& cvec, vec_long& dvec)
1803{
1804   long d_red = (deg(U)+1)/2;
1805
1806   if (IsZero(V) || deg(V) <= deg(U) - d_red) {
1807      return;
1808   }
1809
1810   long du = deg(U);
1811
1812
1813   long d1 = (d_red + 1)/2;
1814   if (d1 < 1) d1 = 1;
1815   if (d1 >= d_red) d1 = d_red - 1;
1816
1817   zz_pXMatrix M1;
1818
1819   ResHalfGCD(M1, U, V, d1, cvec, dvec);
1820   mul(U, V, M1);
1821
1822   long d2 = deg(V) - du + d_red;
1823
1824   if (IsZero(V) || d2 <= 0) {
1825      return;
1826   }
1827
1828   M1(0,0).kill();
1829   M1(0,1).kill();
1830   M1(1,0).kill();
1831   M1(1,1).kill();
1832
1833
1834   zz_pX Q;
1835
1836   append(cvec, LeadCoeff(V));
1837   append(dvec, dvec[dvec.length()-1]-deg(U)+deg(V));
1838   DivRem(Q, U, U, V);
1839   swap(U, V);
1840
1841   ResHalfGCD(M1, U, V, d2, cvec, dvec);
1842
1843   mul(U, V, M1); 
1844}
1845
1846
1847void resultant(zz_p& rres, const zz_pX& u, const zz_pX& v)
1848{
1849   if (deg(u) <= NTL_zz_pX_GCD_CROSSOVER || deg(v) <= NTL_zz_pX_GCD_CROSSOVER) { 
1850      PlainResultant(rres, u, v);
1851      return;
1852   }
1853
1854   zz_pX u1, v1;
1855
1856   u1 = u;
1857   v1 = v;
1858
1859   zz_p res, t;
1860   set(res);
1861
1862   if (deg(u1) == deg(v1)) {
1863      rem(u1, u1, v1);
1864      swap(u1, v1);
1865
1866      if (IsZero(v1)) {
1867         clear(rres);
1868         return;
1869      }
1870
1871      power(t, LeadCoeff(u1), deg(u1) - deg(v1));
1872      mul(res, res, t);
1873      if (deg(u1) & 1)
1874         negate(res, res);
1875   }
1876   else if (deg(u1) < deg(v1)) {
1877      swap(u1, v1);
1878      if (deg(u1) & deg(v1) & 1)
1879         negate(res, res);
1880   }
1881
1882   // deg(u1) > deg(v1) && v1 != 0
1883
1884   vec_zz_p cvec;
1885   vec_long  dvec;
1886
1887   cvec.SetMaxLength(deg(v1)+2);
1888   dvec.SetMaxLength(deg(v1)+2);
1889
1890   append(cvec, LeadCoeff(u1));
1891   append(dvec, deg(u1));
1892
1893
1894   while (deg(u1) > NTL_zz_pX_GCD_CROSSOVER && !IsZero(v1)) { 
1895      ResHalfGCD(u1, v1, cvec, dvec);
1896
1897      if (!IsZero(v1)) {
1898         append(cvec, LeadCoeff(v1));
1899         append(dvec, deg(v1));
1900         rem(u1, u1, v1);
1901         swap(u1, v1);
1902      }
1903   }
1904
1905   if (IsZero(v1) && deg(u1) > 0) {
1906      clear(rres);
1907      return;
1908   }
1909
1910   long i, l;
1911   l = dvec.length();
1912
1913   if (deg(u1) == 0) {
1914      // we went all the way...
1915
1916      for (i = 0; i <= l-3; i++) {
1917         power(t, cvec[i+1], dvec[i]-dvec[i+2]);
1918         mul(res, res, t);
1919         if (dvec[i] & dvec[i+1] & 1)
1920            negate(res, res);
1921      }
1922
1923      power(t, cvec[l-1], dvec[l-2]);
1924      mul(res, res, t);
1925   }
1926   else {
1927      for (i = 0; i <= l-3; i++) {
1928         power(t, cvec[i+1], dvec[i]-dvec[i+2]);
1929         mul(res, res, t);
1930         if (dvec[i] & dvec[i+1] & 1)
1931            negate(res, res);
1932      }
1933
1934      power(t, cvec[l-1], dvec[l-2]-deg(v1));
1935      mul(res, res, t);
1936      if (dvec[l-2] & dvec[l-1] & 1)
1937         negate(res, res);
1938
1939      PlainResultant(t, u1, v1);
1940      mul(res, res, t);
1941   }
1942
1943   rres = res;
1944}
1945
1946void NormMod(zz_p& x, const zz_pX& a, const zz_pX& f)
1947{
1948   if (deg(f) <= 0 || deg(a) >= deg(f)) 
1949      Error("norm: bad args");
1950
1951   if (IsZero(a)) {
1952      clear(x);
1953      return;
1954   }
1955
1956   zz_p t;
1957   resultant(t, f, a);
1958   if (!IsOne(LeadCoeff(f))) {
1959      zz_p t1;
1960      power(t1, LeadCoeff(f), deg(a));
1961      inv(t1, t1);
1962      mul(t, t, t1);
1963   }
1964
1965   x = t;
1966}
1967
1968
1969NTL_END_IMPL
Note: See TracBrowser for help on using the repository browser.