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